Random Forest and Enhanced Sampling Unite: Revealing and Correcting Ghost Errors in Alchemical Free Energy Calculations

随机森林与增强采样联手：揭示并修正炼金术自由能计算中的幽灵误差本文信息标题: 研究炼金术自由能预测中的误差：使用随机森林模型与GaMD 作者: Skanda Sastry and Michael Tae-jong Kim 单位: Genentech Inc, South San Francisco, California, 美国引用格式: Sastry, S., & Kim, M. T.-j. (2025). Investigating Errors in Alchemical Free Energy Predictions Using Random Forest Models and GaMD. Journal of Chemical Information and Modeling. https://doi.org/10.1021/acs.jcim.5c01135 源代码: https://github.com/adnaksskanda/gamdti-paper 摘要当前最先进的抗体-抗原复合物计算结合自由能变化（$\Delta\Delta G$）预测技术，其精度约为$\pm1$ kcal/mol。尽管这对于高通量筛选或亲和力成熟等应用已足够，但对于在临床开发阶段评估翻译后修饰（PTMs）的关键性和影响而言，这一精度仍显不足。那些导致结合能力下降超过50%的PTMs会对实现预期疗效构成重大风险，因此必须严格控制其含量以确保产品质量。50%的解离常数（$K_D$）损失对应于$+0.5$ kcal/mol的$\Delta\Delta G$变化，这意味着计算预测的精度必须达到$\pm0.5$ kcal/mol的阈值，才能在临床阶段具有实际应用价值。在本文中，我们使用常规分子动力学热力学积分（CMD-TI）方法生成$\Delta\Delta G$预测值，并开发了一种结合随机森林（RF）模型和末端态高斯加速分子动力学（GaMD）的误差分析方法。该方法仅需cMD-TI和末端态GaMD数据，即可无偏见地洞察关键自由度（DOF）的采样不足问题。我们发现，大体积侧链的采样不足和关键原子间相互作用的破坏是主要的误差来源，通过我们基于GaMD的误差校正，在误差最大的案例中，预测精度提升超过了1 kcal/mol。当应用于一个包含13个突变的测试集时，基于GaMD的误差校正将均方根误差（RMSE）从$1.06 \pm 0.22$ kcal/mol降低至$0.70 \pm 0.18$ kcal/mol。这项工作不仅开创了利用炼金术自由能预测来评估PTM对生物活性影响的应用，也深入探究了限制其在临床开发中实际应用的关键误差来源。一句话：跑一段GaMD来识别关键DOF，进而指导TI的采样能减小误差。背景治疗性抗体是现代生物医药的基石，其通过与特定抗原的高亲和力结合来发挥治疗作用。在抗体药物的规模化生产过程中，蛋白质不可避免地会发生各种化学修饰，即翻译后修饰（PTMs），如色氨酸氧化、天冬氨酸异构化等。这些PTMs如果发生在抗体-抗原结合界面附近，可能会显著改变结合亲和力，从而影响药物的疗效、药代动力学甚至安全性。因此，准确评估PTMs的影响，并将其作为关键质量属性（CQA）进行严格控制，是生物制药开发中的核心环节。传统上，评估PTM影响主要依赖实验方法，如富集含有特定PTM的抗体亚型，再通过SPR等技术测定其结合活性。然而，这一过程不仅耗时耗力，而且当多种PTMs同时出现时，几乎无法剥离出单一修饰的影响。相比之下，计算模拟方法，特别是炼金术自由能计算，为评估这些点突变或化学修饰对结合自由能的影响（$\Delta\Delta G$）提供了一个高效、精准的理论框架。然而，尽管炼金术自由能计算（如热力学积分TI或自由能微扰FEP）是当前预测相对结合自由能（RBFE）的“金标准”，但其精度仍然存在瓶颈。目前，对于蛋白质-蛋白质相互作用体系，该方法的最佳精度约为$\pm1$ kcal/mol。这一精度水平足以用于抗体亲和力改造的初步筛选，但对于临床阶段的CQA评估，则显得力不从心。一个对产品质量构成严重风险的PTM，其bioactivity影响阈值通常设定为50%，这在热力学上相当于仅仅$+0.5$ kcal/mol的$\Delta\Delta G$变化。因此，计算方法必须达到远超当前水平的$\pm0.5$ kcal/mol精度，才能为临床决策提供可靠依据。这一巨大的“精度鸿沟”是当前领域面临的核心挑战，其背后的误差来源——无论是力场不准、构象采样不足还是计算方案本身的缺陷——亟待被系统性地揭示和解决。 50%解离常数损失意味着什么？在临床上，如果一个PTM导致抗体的生物活性（通常与结合亲和力相关）损失超过50%，则被认为具有高风险。在热力学层面，这意味着结合变得更弱，解离常数$K_D$增大。具体来说，“50%的活性损失”通常指突变体的$K_D$值变为野生型的两倍，即$K_{D,mutant} / K_{D,wildtype} = 2$。根据公式 $\Delta\Delta G = RT \ln(K_{D,mutant} / K_{D,wildtype})$ 在室温下（约298K），这对应于$\Delta\Delta G \approx +0.41$ kcal/mol的变化。为了能够可靠地识别这一变化，计算方法的精度必须显著优于这个值，因此作者提出了$\pm0.5$ kcal/mol的目标。关键科学问题本文旨在解决的核心科学问题是：如何系统性地识别并校正炼金术自由能计算中的微观分子层面采样误差，从而将其预测精度提升至临床应用所需的$\pm0.5$ kcal/mol阈值以下？这不仅仅是一个提升数值精度的问题，更是要深入理解在非物理的炼金术路径中，哪些关键的分子动态行为被错误地表征，并开发出能够“对症下药”的诊断和修正策略。创新点创新的误差诊断框架：首次提出了一种无偏见的（untargeted）误差诊断新方法，该方法巧妙地将机器学习（随机森林）与增强采样（GaMD）相结合，能够从复杂的动力学数据中自动识别出导致计算误差的关键分子自由度（DOF）。揭示核心误差来源：通过该框架，系统性地 pinpoint 了炼金术计算中两个主要的误差来源：一是大体积氨基酸侧链（如Trp）的旋转异构态采样不足；二是在炼金术中间态，由于混合势场的人为效应导致的关键盐桥等原子间相互作用的破坏。精准的误差校正策略：针对上述误差来源，开发了相应的校正方法（如基于GaMD构象分布对TI数据进行过滤或重加权，以及使用距离限制来强制维持关键相互作用），在误差最大的案例中实现了超过1 kcal/mol的精度提升。方法和体系作者采用了一套结合常规MD、增强采样MD和机器学习的综合性方法流程，详见图2。 1. 模拟体系与数据集准备实验数据集：本文使用的基准数据集来源于已发表的文献，主要包括hu4D5-5、mab1和mab2三个抗体系统的一系列单点突变及其对应的实验测定结合能数据。hu4D5-5是人源化抗p185HER2抗体4D5的一个变体，与乳腺癌靶点Erbb2抗原结合。结构准备：抗体-抗原复合物的初始结构来源于PDB数据库（如hu4D5-8的冷冻电镜结构，PDB ID: 6OGE）。hu4D5-5的结构是通过在hu4D5-8上引入两个点突变（VH-V102Y 和 VL-E55Y）构建的。为了节省计算资源，模拟中对抗原蛋白进行了截断，仅保留了靠近结合界面的135个残基。 MD模拟设置：力场与溶剂：所有模拟均采用AMBER20软件包，力场为ff14SB，水模型为TIP3P。体系被溶于一个半径为10 Å的水盒子中，并加入0.15 M的NaCl以模拟生理盐浓度。拓扑构建：使用AmberTools20中的tLEaP和parmed工具准备拓扑文件。对于非天然氨基酸（甲硫氨酸亚砜），使用Gaussian 09和antechamber进行力场参数化。 cMD-TI协议：每个突变计算包含5个重复。体系首先在$\lambda=0.5$下进行能量最小化和升温弛豫，然后进行双向串行平衡，最后在12个$\lambda$窗口下分别进行5 ns的production模拟。每个$\lambda$窗口用于分析的帧数（构象数）为 200帧。 GaMD协议：为了获得更可靠的构象分布，对每个突变的端点态（野生型和突变型）进行了5次重复的、每次300 ns的GaMD增强采样模拟。 2. 随机森林(RF)关键自由度筛选这是本文的核心创新，目的是从海量构象信息中找出导致误差的“罪魁祸首”。详见文末附录。数据集的每一行代表TI模拟过程中的一个单一快照（即一个构象）。对于同一帧，计算机会记录其对应的能量导数值$dV/d\lambda$。特征(Feature)提取：首先，通过GaMD轨迹确定体系的最低能构象簇。然后，在突变位点周围5 Å的球形区域内，定义一系列几何参数作为候选特征，主要包括侧链的二面角（rotamers）和原子间的距离（interatomic distances）。目标变量(Target)定义：RF模型要预测的目标不是原始的能量导数$dV/d\lambda$，而是经过高斯求积权重$w_j$加权后的值，即$w_j \cdot dV/d\lambda$。这使得模型能更直接地关注对最终$\Delta G$积分贡献最大的项。特征筛选与模型训练：使用scikit-learn库进行建模。首先剔除相关性过高（Pearson $r > 0.5$）的冗余特征。然后使用递归特征消除（Recursive Feature Elimination）方法进一步筛选，保留最重要的75%特征。最后，使用这些筛选后的特征训练一个随机森林回归模型，并通过贝叶斯超参数调优来优化模型性能。关键自由度(DOF)识别：模型训练完成后，利用随机森林内置的“基于不纯度的平均特征重要性（mean impurity-based feature importance）”指标，量化每个DOF对预测$w_j \cdot dV/d\lambda$的贡献度。得分最高的DOF即被认为是影响能量计算的关键自由度。 3. 使用的软件工具总结 MD模拟: AMBER20, AmberTools20 (tLEaP, parmed) 增强采样: GaMD 量子化学计算: Gaussian 09 机器学习: scikit-learn 轨迹分析: CPPTRAJ, PyReweighting 分子可视化: VMD 研究内容与结果初始TI预测的性能基准作者首先在一个包含20个有定量实验数据的抗体突变数据集上，评估了他们标准cMD-TI流程的性能。图1：经验ΔΔG与预测ΔΔG的对比图。该图展示了包含所有定量实验结果的案例中，初始TI预测值（纵轴）与实验测量值（横轴）的比较。理想情况下，所有数据点应落在对角虚线上。虽然整体趋势良好（斜率0.788），但均方根误差（RMSE）为0.94 kcal/mol，且许多数据点落在了$\pm1$ kcal/mol的误差区间（点线之间）之外。分析发现，涉及大体积侧链（如Phe, Tyr, Trp）或电荷变化的突变，误差往往更大。创新的RF+GaMD联合误差诊断流程为了剖析这些误差的根源，作者设计了一套创新的诊断流程。图2：TI计算与误差模式分析方法的图形化示意图。该图展示了整个工作流程：(左上) 首先通过常规的TI计算获得初始的$\Delta\Delta G$；(中上) 在突变位点周围5Å的局部环境中测量各种DOF；(右上) 将这些DOF作为输入，加权的$dV/d\lambda$作为输出，训练一个随机森林模型，以识别出对能量影响最大的关键DOF；(中下) 利用GaMD增强采样的轨迹生成这些关键DOF的自由能分布图（PMF）；(左下) 将常规TI模拟对关键DOF的采样情况与GaMD的PMF进行对比，找出采样不一致的地方，并据此进行校正。通过该流程，作者识别出了导致TI计算不准确的关键DOF。跑GaMD不需要事先知道关键DOF？在这个工作流程中，跑GaMD时不需要事先知道哪个或哪些DOF是关键的。这正是该方法“无偏见”（untargeted）的核心优势所在。 GaMD的角色是作为一个独立的、更可靠的“黄金标准”来使用。它通过施加一个偏置势能，对体系的整个势能形貌进行增强采样，目的是尽可能地探索所有可能的构象，并生成一个接近真实平衡态的自由能分布图（PMF）。这个过程是全局性的，不针对任何特定的DOF。关键DOF的识别是在之后发生的。流程是：并行计算：独立地运行常规TI模拟和GaMD增强采样模拟。事后诊断：利用随机森林模型，分析TI轨迹和能量数据，从事后诸葛亮的角度找出哪些DOF对能量计算影响最大。交叉验证：将RF模型找出的关键DOF在TI模拟中的表现，与GaMD这个“黄金标准”进行对比，从而确认采样错误。怎么根据PMF校正采样的？详见附录。 graph TD A["发现TI采样与GaMD PMF不一致"]; A --> B{{"误差类型是什么？"}}; B -- "构象态采样比例错误 (例如：大体积侧链)" --> C1; B -- "关键相互作用持续性破坏 (例如：盐桥断裂)" --> D1; subgraph "方法二：施加限制并重算" direction LR D1["1.从GaMD PMF中 确定关键相互作用的 正常几何范围 (如距离<5Å)"] --> D2["2.根据该范围 设置一个NMR式的距离限制"]; D2 --> D3["3.完全重新进行TI模拟 在所有λ窗口中施加该距离限制"]; D3 --> D4["4.新模拟的结果 即为校正后的ΔΔG"]; end subgraph "方法一：过滤与重加权" direction LR C1["1.从GaMD PMF中 识别有效的低能构象态 (A, B...)" ] --> C2["2.过滤TI轨迹 只保留属于有效构象态的帧"]; C2 --> C3["3.将保留的帧分组 并为每个构象态(A, B...) 单独计算ΔG (ΔG_A, ΔG_B...)"]; C4["4.从GaMD PMF中 获取各构象态的布居比例 (Area_A, Area_B...)"]; C3 & C4 --> C5["5.线性组合得到校正结果 ΔG_corr = Area_A·ΔG_A + Area_B·ΔG_B"]; end 表1：由随机森林模型识别出的误差最大案例中的前5个最重要自由度 rank hu4D5-5 VH-R50A (charging step) hu4D5-5 VH-W95A mab2 VL-Y→R (charging step) mab2 VH-T→Y hu4D5-5 VL-F53N 1 Ag-E71:VH-R50 salt bridge dist VH-W95 chi1 Ag-D161:VL-R49 salt bridge dist Ag-V117:VH-Y53 H-bond dist Ag-C117:VL-N53 H-bond dist 2 VH-R50 chi1 VH-W95 chi2 VL-R49:VL-S50 H-bond dist VL-Y53 chi1 VL-N53 chi1 3 VH-R94 chi4 VH-V48 chi1 VL-S53 chi1 VL-T53 chi1 Ag-M102 chi1 4 VH-F100 chi1 VL-T94:VH-R50 H-bond dist VL-S50 chi1 VL-T53:VL-N51 H-bond dist Ag-N53 chi2 5 Ag-E71 chi3 Ag-E71 chi3 Ag-R157 chi1 VL-Y53 chi2 Ag-N120 chi1 注：表格内容根据原文Table 1整理。加粗的特征是作者后续使用GaMD自由能图进行深入检查的特征。案例分析：揭示三大核心误差来源案例1：大体积侧链采样不足 (Bulky Side-Chain Undersampling) 在hu4D5-5 VH-W95A（色氨酸突变为丙氨酸）这个误差高达1.88 kcal/mol的案例中，RF模型指出，W95侧链的两个二面角（chi1/chi2）是影响能量计算的最关键DOF。图3：(A, C) 完整的和 (B, D) 校正后的TI采样与VH-W95 chi1/chi2旋转角空间的GaMD自由能形貌图的比较，分别对应结合态(A, B)和非结合态(C, D)。图中，背景的彩色热图代表由GaMD增强采样得到的“真实”自由能地貌，其中颜色越深的区域能量越低，是侧链最应该停留的构象。而灰色的散点则代表在常规TI模拟中，侧链实际访问过的构象。在(A)和(C)中可以看到，大量的TI采样点（灰色散点）散落在高能量区域，并未准确地集中在GaMD发现的两个主要低能区域（能量阱）。更重要的是，TI模拟对这两个能量阱的采样比例（例如在结合态，TI采样比例为33.6% vs 66.4%）与GaMD计算的真实布居比例（GaMD: 24.5% vs 75.5%）存在显著偏差。作者通过过滤掉无效的TI采样帧，并根据GaMD的比例对两个能量阱的贡献进行重新加权，最终将预测误差从1.88 kcal/mol降至0.44 kcal/mol。案例2：关键盐桥相互作用的破坏 (Violation of Key Salt Bridge Interactions) 在hu4D5-5 VH-R50A和mab2 VL-Y→R这两个涉及电荷变化的突变案例中，RF模型发现，影响计算的最关键DOF是抗体与抗原之间的一个关键盐桥的距离。图5：(A) hu4D5-5复合物中的Ag-E71:VH-R50A盐桥和(B) mab2复合物中的Ag-D161:VL-Y→R盐桥的可视化。图中展示了由RF模型识别出的关键盐桥。图4：(A-C) hu4D5-5 VH-R50A和(D-F) mab2 VL-Y→R中关键盐桥距离的GaMD自由能形貌图（曲线）与TI采样（直方图）的对比。(A, D)为原始TI模拟，(B, E)为过滤后的TI数据，(C, F)为施加NMR距离限制后的TI模拟。在(A)和(D)中，可以惊人地发现，在大部分TI模拟帧中（绿色直方图），该盐桥的距离都远远超过了5Å，说明这个关键的相互作用在计算过程中被人为地破坏了。作者推测这是由于炼金术中间态的混合势场削弱了静电作用所致。通过(B, E)过滤掉盐桥破坏的帧，或(C, F)在TI模拟中施加距离限制来强制维持盐桥，预测精度都得到了超过1 kcal/mol的显著提升。小编补充：从图4来看，似乎过滤数据分布也差不多，但很可能普通TI散掉了就一直散掉了，采不了多少数据？还是约束着好。案例3：关键氢键相互作用的破坏 (Violation of Key Hydrogen Bonds) 对于mab2 VH-T→Y和hu4D5-5 VL-F53N这两个案例，RF模型识别出的关键DOF是分子间的氢键距离。与前两种情况类似，TI模拟也未能准确捕捉这些氢键的正确构象。然而，对这类更动态、更复杂的相互作用进行校正要困难得多。例如，在mab2 VH-T→Y中，由于自由能形貌呈现双峰分布，简单的加权求和难以实现。在hu4D5-5 VL-F53N中，施加距离限制甚至导致了体系无法正确平衡。这表明，虽然该框架能有效识别问题，但对某些复杂情况的修复仍是未来的挑战。总体校正效果表2：对误差最大的几个预测进行GaMD校正的效果总结 | case | original $\Delta\Delta G$ (kcal/mol) | corrected $\Delta\Delta G$ (kcal/mol) | empirical $\Delta\Delta G$ (kcal/mol) | original error (kcal/mol) | corrected error (kcal/mol) | model R² | theorized source of error | | :— | :— | :— | :— | :— | :— | :— | :— | | hu4D5-5 W95A | $3.19 \pm 0.71$ | $4.63 \pm 0.47$ | 5.07 | 1.88 | 0.44 | 0.53 | bulky side-chain undersampling | | hu4D5-5 R50A | $2.98 \pm 1.02$ | $5.14 \pm 1.05$ | 4.58 | 1.60 | 0.56 | 0.73 | salt bridge violation | | mab2 VL-Y→R | $0.61 \pm 1.04$ | $-0.87 \pm 1.01$ | < -0.83 | > 1.43 | - | 0.48 | salt bridge violation | | mab2 VH-T→Y | $1.69 \pm 0.75$ | $0.61 \pm 0.73$ | 0 | 1.69 | 0.61 | 0.76 | hydrogen bond violation | | hu4D5-5 VL-F53N| $-0.56 \pm 0.57$ | - | 1.19 | 1.75 | - | 0.64 | hydrogen bond violation | 注：表格内容根据原文Table 2整理。不确定度为95%置信区间。original：普通TI；empirical：实验世界最终，作者将此校正方法应用到包含13个突变的整个hu4D5-5数据集，进一步验证了其普适性。表3：对hu4D5-5数据集($n=13$)的误差分析：原始方案、延长模拟方案与RF+GaMD校正方案的对比 | protocol | MAE | RMSE | | :— | :— | :— | | original protocol (5 ns per $\lambda$) | $0.82 \pm 0.18$ | $1.06 \pm 0.22$ | | extended protocol (25 ns per $\lambda$) | $0.71 \pm 0.18$ | $0.93 \pm 0.23$ | | RF + GaMD correction | $0.53 \pm 0.16$ | $0.70 \pm 0.18$ | 注：表格内容根据原文Table 3整理。MAE为平均绝对误差，RMSE为均方根误差。结果表明，简单地将模拟时间延长5倍，对精度的提升有限，而RF+GaMD校正方法则取得了实质性的改进，使RMSE非常接近$\pm0.5$ kcal/mol的目标精度。深入讨论本文的讨论部分对研究的发现及其意义进行了深刻的阐述。盐桥破坏是炼金术计算中的一个普遍且棘手的问题：作者强调，在炼金术中间态，混合势场会系统性地削弱静电相互作用，导致关键盐桥的“局部解离”。这是一个方法本身的缺陷，无法通过简单延长模拟时间来解决。作者将他们基于物理的距离限制校正方法与文献中其他经验性的校正方案进行对比，认为他们的方法虽然计算成本更高，但更为严谨和可靠，特别适用于对精度要求极高的场景。机器学习赋能无偏见的误差诊断：本文最大的方法学创新在于利用RF模型实现了一种“无偏见”的误差溯源。在不具备任何先验知识的情况下，该模型能自动从纷繁的动力学数据中识别出对能量计算起决定性作用的少数几个DOF。这为解决自由能计算中的“未知之不知”问题提供了一个强大的工具。对已知问题的再验证：RF模型能够自动识别出“大体积侧链采样不足”这一领域内公认的难题，这本身就强有力地验证了该诊断框架的有效性。作者指出，他们的框架可以作为更高级的增强采样方法（如ACES）的前导步骤，为其指明需要增强采样的关键DOF，从而提高效率。对力场误差的评估：在经过一系列采样校正后，预测的RMSE降低到了0.70 kcal/mol。考虑到实验测量本身也存在约$\pm0.3$ kcal/mol的误差，这意味着由ff14SB力场本身带来的误差可能非常小（约0.2-0.4 kcal/mol）。这表明，在当前阶段，改善采样问题比优化力场参数对提升预测精度的贡献可能更大。 Q\&A Q1: 为什么简单地延长TI模拟时间通常无法修复这些预测误差？ A1: 因为这些误差很多是系统性误差（systematic error），而非随机误差。例如，在电荷变化突变中，炼金术路径中间态的“混合势场”会人为地削弱静电相互作用。这导致关键的盐桥即使在理论上应该存在的情况下也容易断裂。无论模拟时间多长，只要这个势场本身存在缺陷，盐桥就可能一直处于被破坏的状态。这并非采样不足的问题，而是方法本身的“人造缺陷”（artifact）。 Q2: 随机森林模型（RF）在其中扮演了什么角色？为什么不直接分析所有可能的分子自由度（DOF）？ A2: RF模型扮演了“筛子”或“侦探”的角色。在一个复杂的蛋白质体系中，分子自由度（如所有侧链的旋转角、所有原子间的距离）的数量是巨大的。绝大多数DOF的变化对我们关心的$\Delta\Delta G$计算影响甚微，它们是“噪音”。RF模型通过监督学习，能够从海量的DOF中，找出与能量导数$dV/d\lambda$相关性最强、即对最终结果影响最大的那几个“关键自由度”。这使得后续的分析可以集中火力解决主要矛盾，而不是在无关的噪音中大海捞针。 Q3: GaMD增强采样给出的自由能分布就一定是“正确”的吗？这个方法的核心假设是什么？ A3: GaMD并不保证绝对“正确”，但它是一种增强采样方法，相比于几纳秒的常规MD（cMD），它能更快速、更广泛地探索分子的构象空间，因此其得到的自由能分布更有可能接近体系的真实平衡态分布。该方法的核心假设是：由长时间GaMD模拟得到的构象能量地貌，比短时间的常规TI模拟所采样的构象，能更准确地反映体系的真实热力学性质。当然，使用GaMD时，研究者也需要通过检查收敛性等方式来确认其结果的可靠性。 Q4: 论文中使用距离限制（restraint）来修复盐桥断裂问题，但施加限制本身不会引入新的能量项，从而影响自由能计算吗？ A4: 问得非常好，这确实是一个严谨性问题。是的，施加限制会改变系统的哈密顿量，理论上需要计算并扣除这个限制所贡献的自由能。在一些体系中（如小分子-蛋白），有成熟的方法（如Boresch restraints）来解析地计算这一项。但在复杂的蛋白-蛋白界面，这个问题尚无标准解法。作者在文中也承认了这一点，他们认为，由盐桥断裂引入的巨大误差（>1 kcal/mol）远大于忽略限制自由能所带来的微小误差，因此在当前阶段，这是一个实用且有效的近似处理方法。 Q5: 这个RF+GaMD框架与其他改进炼金术计算的方法（如ACES）有何不同？ A5: 它们是互补而非竞争的关系。像ACES（炼金术增强采样）这类方法，旨在通过在炼金术路径上对某些“慢”自由度进行增强采样来加速收敛。但一个前提是，你需要预先知道应该对哪些自由度进行增强采样。而本文提出的RF+GaMD框架的核心贡献之一，就是提供了一种无偏见的、自动化的方法来识别出这些需要被特别关注的关键自由度。因此，可以将该框架视为ACES等更高级采样方法的前置步骤：先用RF+GaMD做“侦查”，找出问题所在，再用ACES等方法进行“精确打击”。关键结论与批判性总结关键结论本文提出并验证了一个结合随机森林（RF）和高斯加速分子动力学（GaMD）的创新框架，该框架能够以一种无偏见（untargeted）的方式，系统性地诊断和校正炼金术自由能计算中的采样误差。研究发现，大体积侧链（如色氨酸）的旋转异构态采样不足和在炼金术中间态由于混合势场导致的关键盐桥人为断裂，是导致$\Delta\Delta G$预测不准确的两大主要来源。通过应用基于GaMD的针对性校正策略（对不同构象态进行重加权或施加距离限制），在误差最大的案例中，预测精度提升超过1 kcal/mol。对于包含13个突变的hu4D5-5数据集，该方法将整体均方根误差（RMSE）从$1.06 \pm 0.22$ kcal/mol降至$0.70 \pm 0.18$ kcal/mol，非常接近临床应用所需的0.5 kcal/mol精度目标。研究明确表明，简单地将模拟时间延长5倍（从每个$\lambda$窗口5 ns增加到25 ns）对精度的提升非常有限（RMSE仅从1.06 kcal/mol降至0.93 kcal/mol），这证实了误差主要来源于系统性缺陷而非随机采样不足，因此必须采用更具针对性的校正方法。批判性总结潜在影响：该工作为突破炼金术自由能计算的精度瓶颈提供了一个强大、系统且思路清晰的框架。其“无偏见的误差识别能力”尤为重要，能够在使用常规TI计算的基础上，为更高级的增强采样方法（如ACES）指明需要关注的关键自由度。这有望将高精度自由能计算从少数专家的“炼丹术”推广为更可靠、更自动化的标准流程，从而加速其在治疗性抗体药物临床前CQA评估等要求苛刻的工业领域的应用。存在的局限性：尽管对盐桥断裂的校正效果显著，但对于更瞬态、更复杂的相互作用（如氢键网络）的校正仍具挑战性，因为这些相互作用的自由能形貌可能呈现复杂的多峰分布，难以直接修复。此外，施加距离限制所贡献的自由能未能被严格计算，这在理论上仍是一个待解决的问题。同时，该框架无疑增加了额外的计算成本，因为它需要在标准TI计算之外进行长时间的GaMD模拟和机器学习模型训练。未来研究方向：未来的工作可以致力于将该框架推广到更多的PTMs类型（如天冬酰胺脱氨等）和更复杂的分子体系中。同时，开发更普适、更自动化的方法来校正氢键等瞬态相互作用，以及探索如何在蛋白-蛋白体系中严格处理限制自由能的贡献，将是该领域的重要发展方向。附录1：随机森林(RF)关键自由度筛选高斯求积权重$w_j$的来源：高斯求积是一种经典的数值积分方法，其目的是用一个离散的加权和来高度精确地近似一个连续积分。对于热力学积分中自由能的计算，其基本形式为： $\Delta G = \int_{0}^{1} \left\langle \frac{\partial V(\lambda)}{\partial \lambda} \right\rangle_{\lambda} d\lambda$ 为了在计算机中求解，该积分被离散化。在使用N点高斯求积法时，这个积分被近似为在N个特定的$\lambda$坐标点（$\lambda_j$）上计算值的加权和： $\Delta G \approx \sum_{j=1}^{N} w_j \cdot \left\langle \frac{\partial V(\lambda)}{\partial \lambda} \right\rangle_{\lambda_j}$ 在这项工作中，作者使用了12点高斯求积方案（即$N=12$）。这些$\lambda_j$坐标点和它们对应的权重$w_j$是根据高斯求积的数学理论预先确定的，旨在最大化数值积分的精度。该论文将这些权重作为标准数值方法的组成部分直接使用，并未详细阐述其数学推导过程。递归特征消除（Recursive Feature Elimination, RFE）的具体操作：该方法通过一个迭代过程来系统性地减少特征数量，以找到性能最佳的特征子集。其操作流程如下：初始训练：使用全部$p$个候选特征$S_0 = {f_1, f_2, \dots, f_p}$训练一个基础模型（本文中是一个最大深度为5的决策树回归器）。重要性排序：根据训练好的模型，计算并排序所有特征的重要性。特征剔除：移除最不重要的预设百分比的特征。根据论文描述，每次迭代剔除5%的特征。这个过程可以表示为： $S_{k+1} = S_k \setminus R_k$ 其中，$S_k$是第$k$次迭代的特征集，$R_k$是该轮中被评定为最不重要的5%特征的集合。循环迭代：重复步骤1至3，直到特征数量达到预设的目标。在本文中，该过程持续进行，直到保留原始特征集中最重要的75%为止。贝叶斯超参数调优的具体操作：这是一种用于自动寻找机器学习模型最佳超参数组合的优化算法。其目标是找到一组能使模型性能最大化的超参数$\theta^*$。定义目标函数：首先定义一个评估模型性能的目标函数$J(\theta)$。在本文中，目标函数被设定为5折交叉验证后的平均$R^2$值，这代表了模型的泛化能力。构建代理模型：贝叶斯优化使用一个概率模型（通常是高斯过程）来拟合已观察到的超参数点$(\theta, J(\theta))$与目标函数之间的关系。优化搜索：整个搜索过程共进行200次迭代。前50次通过拉丁超立方采样进行随机探索，以获得对超参数空间的初步了解。随后的150次迭代则由贝叶斯模型指导，通过一个“采集函数”来智能地选择下一个最有希望提升性能的超参数组合进行尝试，从而高效地找到全局最优解。整个优化问题可表示为： $\theta^* = \arg\max_{\theta \in \Theta} J(\theta)$ 其中$\Theta$是所有可能的超参数组合空间。基于不纯度的平均特征重要性的具体计算：这是决策树和随机森林模型中常用的一种评估特征重要性的方法。对于回归任务，其核心是计算每个特征对“方差减少”的贡献。节点方差：对于树中的任意一个节点$m$，其包含的数据点的方差定义为： $\text{Var}(m) = \frac{1}{N_m} \sum_{i \in \text{node } m} (y_i - \bar{y}_m)^2$ 其中$N_m$是节点$m$中的样本数，$y_i$是样本值，$\bar{y}_m$是节点内所有样本的平均值。分裂带来的方差减少：如果一个节点$m$使用特征$f$进行分裂，产生左右两个子节点，那么这次分裂带来的方差减少量（即该节点的重要性）为： $\Delta I(m, f) = \text{Var}(m) - \left( \frac{N_{\text{left}}}{N_m} \text{Var}(\text{left}) + \frac{N_{\text{right}}}{N_m} \text{Var}(\text{right}) \right)$ 特征在单棵树中的重要性：一个特征$f$在单棵决策树$T$中的总重要性，是它在所有用它进行分裂的节点上带来的方差减少量的总和。特征在森林中的重要性：在随机森林中，一个特征的最终重要性是它在所有树中的重要性的平均值。为了结果的稳健性，作者通过5次重复的5折交叉验证共训练了25个模型，最终的特征重要性是这25个模型计算出的重要性分数的平均值。附录2：校正采样的细节怎么根据PMF校正采样的？根据识别出的不同误差来源，作者采用了两种不同的、具有针对性的校正策略： 1. 针对构象态采样比例错误的校正（过滤与重加权）这种方法主要用于处理像大体积侧链采样不足（如W95A案例）这样的问题，即TI模拟虽然找到了正确的低能构象态，但对它们的采样比例是错误的。第一步：识别构象态。首先，根据GaMD生成的PMF，确定体系存在几个主要的低能构象微观态（microstates）以及它们各自的能量盆。例如，在W95A案例中，GaMD发现W95侧链主要存在两个稳定的旋转异构态。第二步：过滤TI数据。检查常规TI模拟的每一帧，将所有不属于GaMD识别出的任何一个低能构象态的帧全部过滤掉、丢弃。这些被认为是物理意义不大或采样错误的“噪音”数据。第三步：分别计算各态的自由能。对于过滤后剩下的数据，将其按照所属的构象态进行分组。然后，为每一个构象态单独计算其炼金术自由能变化$\Delta G$。例如，如果存在两个微观态，就会得到$\Delta G_1和\Delta G_2$。第四步：根据GaMD比例进行重加权。最后，根据GaMD的势能面（Potential of Mean Force, PMF）计算出各个微观态的真实热力学布居比例（即自由能盆的面积或体积占比，例如$\text{Area}_1$和$\text{Area}_2$）。用这个比例作为权重，对上一步分别计算出的自由能进行线性组合，得到最终校正后的总自由能： $\Delta G_{\text{corrected}} = (\text{Area}_1 \times \Delta G_1) + (\text{Area}_2 \times \Delta G_2) + \dots$ 这个过程本质上是用热力学积分（Thermodynamic Integration, TI）的局部能量信息，结合增强采样分子动力学（GaMD）的全局构象分布信息，来重构一个更准确的自由能值。 2. 针对关键相互作用破坏的校正（施加距离限制）这种方法主要用于处理像关键盐桥断裂（如R50A案例）这样的问题，即TI模拟系统性地无法采样到某个本应存在的关键相互作用。第一步：识别相互作用。通过GaMD的PMF确认某个关键相互作用（如盐桥）在平衡态下是稳定存在的，并确定其正常的距离范围（例如< 5 Å）。第二步：施加距离限制并重新模拟。作者发现，简单地过滤数据会导致样本量急剧下降。因此，他们采用了一种更稳健的方法：重新进行一次TI模拟。在这次新的模拟中，他们施加了一个NMR式的距离限制（distance restraint），强制性地将形成盐桥的两个原子基团的距离约束在GaMD确定的合理范围内。第三步：使用限制性模拟的结果。这个限制有效地阻止了盐桥在炼金术中间态的人为断裂，确保了这一关键相互作用在整个计算过程中的完整性。最终的$\Delta\Delta G$直接采用这次限制性TI模拟的结果。虽然从理论上讲，施加限制本身会引入额外的自由能项，但作者认为，由盐桥破坏引入的巨大误差（>1 kcal/mol）远大于忽略限制自由能所带来的微小误差，因此这是一个在实践中非常有效的校正策略。如何从GaMD PMF中确定关键相互作用的正常几何范围？从GaMD（高斯加速分子动力学）生成的PMF（Potential of Mean Force，平均力势）图中确定相互作用的正常几何范围，主要依赖于对自由能形貌的解读。这个过程可以分为两步：第一步：生成并观察自由能分布图首先，需要针对感兴趣的几何参数（例如形成盐桥的两个原子团之间的距离）运行GaMD模拟并计算其一维PMF。这个PMF图的纵轴是自由能，横轴是距离。一个热力学稳定的相互作用会在图中表现为一个清晰、深刻的能量阱（energy well）。在论文的图4中，这个能量阱体现为相对丰度（Relative Abundance）图上的一个尖锐、高耸的山峰。第二步：根据能量阱定义范围 “正常几何范围”就是这个能量阱所覆盖的距离区间。具体操作是：定位能量最低点：找到能量阱最深处（即概率峰值最高处）对应的距离值。这代表了该相互作用最可能存在的距离。确定边界：从能量最低点向两侧延伸，直到自由能开始急剧上升的位置为止。这个能量急剧上升的“井壁”就定义了稳定相互作用的边界。应用临界值：在实践中，可以根据物理化学常识设置一个合理的临界值（cutoff）。例如，对于盐桥，通常认为带电原子团之间的距离在4-5 Å以内才算形成有效的相互作用。在论文的图4中，GaMD的PMF清晰地显示能量阱完全位于5 Å以内，因此作者采用“距离 < 5 Å”作为判断盐桥是否完整的标准是合理且有数据支持的 2。附录3：SI的信息 1. 完整的实验基准数据集 (Table S1) SI提供了用于验证计算结果的全部23个突变的完整实验数据。这包括每个突变的来源文献、实验方法（如SPR、滴定量热法）、原始测量值（如Kd值），以及最终转换为$\Delta\Delta G$ (kcal/mol)的结果。文件还澄清了数据处理的细节，例如在hu4D5-5数据存在多个报告值时，优先选择SPR数据，但对于解离速率过快的突变（如W95A），则根据与原作者的沟通改用等温滴定微量热法（ITC）的数据。 2. 完整的初始TI计算结果 (Table S3) 与实验数据相对应，SI列出了所有23个突变的初始TI计算预测值（$\Delta\Delta G$）及其不确定度。该表格还对每个突变进行了分类，明确标注了其是否涉及大体积侧链（bulky side chain）、电荷变化（charge-changing）或两者兼有。这使得读者可以直接比较不同类型突变的预测难度和误差大小。 3. 误差来源的排他性证据 (Table S2) 在分析涉及电荷变化的突变时，炼金术转化通常分为范德华（vdW）和静电（charging）两个步骤。主文假设误差主要来源于静电步骤。 Table S2提供了关键的“排除法”证据：当作者将RF+GaMD校正方法仅应用于误差最大的两个电荷变化突变（R50A和Y→R）的范德华步骤时，预测精度的改善微乎其微（trivial change）。这有力地证明了误差几乎完全集中在静电（charging）步骤，与主文中观察到的盐桥破坏现象高度吻合。 4. 随机森林（RF）模型的详细参数与定义 (Table S4, S4.2) 为了提高研究的可复现性，SI提供了RF分析的更多细节。Table S4列出了主文中提到的前5个最重要自由度（DOF）的定量重要性分数。 S4.2节提供了每个关键DOF的精确原子定义。例如，它明确定义了“Ag-E71:VH-R50盐桥距离”是“抗原E71残基的CD原子与抗体VH链R50残基的CZ原子之间的距离”。这些精确的定义对于其他研究者复现或借鉴该方法至关重要。 5. 方法的稳健性验证 (Table S5, Figures S1-S4) 为了排除误差是由于特定的“两步法”电荷转化方案引起的可能性，作者使用了一种更新的“一步法”转化方案（使用smoothstep软核势）重新计算了两个关键的电荷变化突变。结果显示，即使在“一步法”中，同样的盐桥破坏问题依然存在。并且，施加距离限制同样能有效地校正误差。这表明盐桥破坏是一个与炼金术混合势场相关的普遍性问题，而非特定计算方案的产物。 6. 发现的普适性——对外部数据的验证 (Figures S5-S11) 为了验证其发现的普适性，作者将其分析思路应用到了一个完全不同的、已发表的barstar-barnase蛋白复合物体系中，该体系的某些突变在原研究中也存在无法解释的巨大误差。作者对这些出错的突变进行了GaMD模拟，结果发现，在每一个出错的案例中，都存在一个先前未被讨论的关键盐桥或氢键相互作用。这强烈暗示，这些外部数据集中的误差很可能也是由同样的关键相互作用破坏机制导致的，从而极大地增强了本文结论的普适性。 7. 对比“增加算力”与“智能校正”的效果 (Table S6, S7) SI提供了最有说服力的数据之一：简单粗暴地增加算力是否能解决问题？作者将所有模拟的采样时间增加了5倍（从每个λ窗口5 ns延长到25 ns）。结果显示，5倍的算力投入对精度的提升非常有限（RMSE仅从1.06轻微降至0.93 kcal/mol），甚至在某些情况下预测结果反而变得更差。相比之下，应用RF+GaMD智能校正方法的RMSE则显著降低至0.70 kcal/mol。这组对比有力地证明了文中所述的误差是系统性误差，无法通过“大力出奇迹”来解决，必须采用本文提出的这种更智能的诊断和校正策略。

Molecular Dynamics · 2025-08-22

VMD Gets a New Tool! PACKMOL-GUI: One-Stop Solution for Complex Molecular System Building

VMD再添利器！PACKMOL-GUI：一站式搞定复杂分子体系的“搭积木”难题本文信息标题: PACKMOL-GUI: An All-In-One VMD Interface for Efficient Molecular Packing 作者: Jian Huang, Chenchen Wu, Xiner Yang, Zaixing Yang, Shengtang Liu, Gang Yu 单位: Soochow University, Children’s Hospital of Zhejiang University School of Medicine 引用格式: Huang, J., Wu, C., Yang, X., Yang, Z., Liu, S., & Yu, G. (2025). PACKMOL-GUI: An All-In-One VMD Interface for Efficient Molecular Packing. Journal of Chemical Information and Modeling, 65, 778-784. 摘要 PACKMOL是计算化学领域广泛使用的分子建模工具。然而，长期以来，它一直缺乏一个强大的、集参数设置与分子和几何约束可视化于一体的开源图形用户界面（GUI），这在很大程度上阻碍了其巨大优势的发挥。为了解决这一局限，我们开发了一款名为PACKMOL-GUI的VMD插件，它利用了Tcl/Tk工具包的动态可扩展性。该GUI允许用户通过一个直观的面板配置PACKMOL的所有参数，同时借助VMD软件，能够方便地可视化分子结构以及包括立方体、盒子、球体等在内的各种几何约束。VMD与PACKMOL之间的无缝交互，为构建复杂的分子系统提供了一个直观、高效的一体化平台。背景分子动力学（MD）模拟是研究复杂分子系统热力学和动力学行为的核心计算方法。在MD模拟工作流程中，一个至关重要的前提步骤是构建一个包含多种分子混合物的、合理的初始构象。想象一下，要在一个模拟盒子中搭建一个复杂的细胞膜体系，你需要精确地放置成百上千个脂质分子、水分子，甚至还有蛋白质和离子，这就像是在一个微观世界里玩一个极其精密的“搭积木”游戏。为了解决这个分子“堆叠”或“填充”的问题，PACKMOL应运而生，并成为该领域应用最广泛的程序之一。它允许用户在定义的空间区域内（如球体、立方体或更复杂的形状）放置指定数量的不同类型的分子，同时避免原子间的严重重叠。然而，PACKMOL的强大功能长期以来被其原始的命令行操作方式所束缚。用户需要手动编写包含大量坐标、几何约束和分子类型的文本输入文件，这个过程不仅繁琐、耗时，而且极易出错。更重要的是，用户无法直观地看到自己设置的几何约束区域与分子之间的关系，只能在运行结束后通过可视化软件检查结果，这使得调试过程非常低效。尽管之前有研究者尝试开发PACKMOL的GUI，例如GEMS-Pack和Atomistica.online，但它们仍存在诸多不足。GEMS-Pack目前已无法访问，并且其依赖的Python 2.7和PyQt5技术栈面临被淘汰的风险，给安装带来挑战。而Atomistica.online则在PACKMOL参数设置、分子与几何约束的可视化方面功能有限，并且有计算时间限制。因此，科研社区迫切需要一个友好的、开源的、并且能将参数设置、分子可视化和约束可视化三者无缝集成的GUI工具。关键科学问题本文旨在解决的核心科学问题是：如何为功能强大但操作繁琐的PACKMOL程序开发一个稳定、开源且功能全面的一体化图形用户界面，使其能够无缝集成到主流的分子可视化软件（如VMD）中，从而将复杂的命令行输入文件生成过程，转变为一个直观的、“所见即所得”的交互式建模体验，最终大幅提升构建复杂分子体系的效率和便捷性？创新点 VMD插件形式：利用VMD广泛的用户基础及其通过Tcl/Tk脚本的动态可扩展性，将PACKMOL的功能直接集成到科研人员熟悉的可视化环境中，无需修改VMD源码或重新编译。一体化平台：首次实现了一个集参数配置、分子结构可视化和几何约束实时可视化于一体的完整工作流。用户可以直接在VMD窗口中看到设置的几何形状（如球体、盒子），极大地增强了操作的直观性。用户友好设计：提供了丰富的内置功能以提升效率，包括一个包含常用分子（脂质、溶剂、离子等）的共享数据库，以及基于体积或表面积自动估算最大可容纳分子数的功能。开源与跨平台：该工具是开源的，并且由于VMD本身支持Windows、Linux和macOS，PACKMOL-GUI也天然地支持这些主流操作系统。研究内容核心方法：PACKMOL-GUI工作流详解 PACKMOL-GUI的设计遵循PACKMOL程序本身的数据流逻辑，将整个建模过程分解为一系列有序的步骤。用户在VMD的“Extensions”菜单中启动插件后，便可进入其主界面。图1：PACKMOL-GUI工作流概览整个工作流程可以清晰地划分为几个核心模块，从通用参数的初始化开始，到分子导入、空间约束定义，最终生成输入文件并运行PACKMOL。 graph TD direction LR subgraph "PACKMOL-GUI 核心工作流" A("VMD Main Extensions->PACKMOL") --> B("初始化通用参数"); subgraph "通用参数" direction LR C["PACKMOL路径 公差/文件类型/pbc 输出目录等"] end B -- "设置" --> C; B --> D("导入分子"); subgraph "分子数据库" direction LR E[("可用数据集")] end D -- "从数据库加载" --> E; D --> F("设置分子数量"); F --> G("定义空间约束"); subgraph "几何约束可视化" direction LR H["球体/椭球体 圆柱/平面/盒子 高斯曲面"] end G -- "实时显示几何形状" --> H; G --> I("生成输入文件 并运行PACKMOL"); I --> J("输出文件"); end 图2：PACKMOL-GUI的布局 PACKMOL-GUI的界面布局遵循自上而下的逻辑顺序，分为五个核心模块，每个模块由不同颜色的虚线边框明确区分。通用参数模块 (General Parameters Module): 首次使用时，用户需要指定本地PACKMOL程序的可执行文件路径。该模块允许设置全局参数，如公差（tolerance）、输出文件类型（filetype）、周期性边界条件（PBC）等。所有设置（如输出目录、参数等）都会被保存在一个名为packmol_info.json的文件中，方便下次使用。为了方便用户，界面右侧还内嵌了PACKMOL的用户手册，可随时查阅。分子导入模块 (Molecule Import Module): 用户可以通过“Import”, “Delete”, “Refresh”按钮来导入、删除或同步分子列表。该模块集成了一个包含常用生物分子、溶剂、气体分子、离子和纳米材料的数据库，极大地便利了复杂系统的建模。例如，离子类别甚至包括了放射性核素离子。一个关键特性是自动估算最大分子数。我们知道，在一个有限的空间里能塞进多少分子是有限的。PACKMOL-GUI提供了两种估算方法：体积估算法 \[N_{vmax}=\frac{V_{constraints}}{V_{molecule}}\] 公式的通俗解释这个公式用于估算在一个给定的约束体积 $V_{constraints}$ 中，最多可以填充多少个分子。$N_{vmax}$ 是最大分子数，$V_{molecule}$ 是单个分子的体积。这个体积值可以通过MoloVol等工具计算得出。表面积估算法（针对膜系统） \[N_{smax}=\frac{S_{constraints}}{APL_{molecule}}\] 公式的通俗解释对于脂双层这样的膜系统，更关心的是在膜的表面能铺多少个脂质分子。$N_{smax}$ 是最大脂质分子数，$S_{constraints}$ 是约束形状提供的膜表面积，$APL_{molecule}$ 是每个脂质分子的平均占用面积（Area Per Lipid）。约束模块 (Constraints Module): 这是PACKMOL程序最具特色的功能，也是该GUI的核心。用户可以为导入的分子或其中的特定原子添加、修改或删除约束。位置约束: 可以定义分子位于某个几何形状的“内部(inside)”、“外部(outside)”、“上方(over)”或“下方(below)”。几何类型: 支持多种几何形状，包括立方体、盒子、球体、椭球体、平面、圆柱体和高斯曲面。实时可视化: 当用户输入几何参数并按下回车键后，相应的几何形状会立即在VMD的主显示窗口中被绘制出来。用户还可以通过界面上的单选按钮控制形状和标签的显示/隐藏，并修改线条粗细、颜色等，实现了真正的“所见即所得”。输入文件生成与执行模块 (Input File Generation and Execution Module): 在所有参数配置完成后，点击“generate”按钮，即可在左侧的文本框中看到生成的PACKMOL输入文件。用户可以点击“save”保存该文件，同时为了防止文件丢失，程序在生成时会自动在工作目录下保存一个带时间戳的副本。确认无误后，点击“run”按钮即可在后台调用PACKMOL程序执行计算。输出日志模块 (Output Log Module): PACKMOL程序的实时运行状态和输出信息会被重定向到该模块的文本框中，方便用户监控执行过程并快速定位和修正输入文件中的错误。案例研究为了展示PACKMOL-GUI的强大性能，作者复现了两个复杂的分子体系构建任务。案例一：构建双层棕榈酸球形囊泡这是一个来自PACKMOL官网的经典案例，目标是构建一个被水溶液包围的、内部也含有水核的脂质囊泡。图3：内外均有水的双层球形囊泡示例这个复杂的体系需要对水分子和棕榈酸分子施加四种不同的空间几何约束。内部水核 (water-0)：被约束在一个半径为13 Å的球体内部。内层脂质 (palmitoyl-1)：其亲水头部被约束在一个半径为14 Å的球内，而疏水尾部则被约束在一个半径26 Å的球外。外层脂质 (palmitoyl-2)：其疏水尾部被约束在一个半径29 Å的球内，而亲水头部则被约束在一个半径41 Å的球外。外部溶剂 (water-3)：被约束在一个边长为90 Å的立方体盒子内部，同时还要满足位于半径为43 Å的球体外部的条件。在PACKMOL-GUI中，用户可以直观地看到这几个层层相套的球形和立方体约束（如图3a所示），并使用Molcontroller工具将不同分子移动到各自的几何区域内进行预览，从而确保约束设置的准确性。案例二：阳离子MOF材料富集放射性离子这个案例来自作者之前的研究，目标是构建一个包含阳离子金属有机框架（MOF）材料SCU-103、多种竞争性阴离子（OH⁻, NO₃⁻, SO₄²⁻, ⁹⁹TcO₄⁻）、抗衡离子和大量水分子的复杂体系。作者提到，在之前的工作中，他们使用GROMACS和Molcontroller等工具迭代构建这个体系，过程非常繁琐耗时。图4：用于吸附⁹⁹TcO₄⁻的阳离子MOF SUC-103 使用PACKMOL-GUI，这个过程变得异常高效。 MOF约束：首先将SCU-103材料放置在由一个蓝色盒子定义的中心区域。离子约束：在MOF表面的上下两侧，使用黄色和橙色的盒子来定义各种离子的初始分布区域。溶剂约束：最后，使用一个赭石色的盒子来定义整个水溶剂的边界。通过GUI的可视化功能，用户可以清晰地看到代表不同约束区域的彩色盒子（如图4a所示），从而快速、准确地完成整个复杂系统的初始构象搭建。 Q\&A Q1: PACKMOL-GUI相比于之前的GEMS-Pack等GUI工具有哪些本质上的优势？ A1: 最核心的优势是深度集成与可视化。PACKMOL-GUI是作为VMD的插件运行的，这意味着它能直接利用VMD强大的分子可视化和操作能力。用户在设置几何约束时，可以实时在VMD窗口中看到这些约束（如球体、盒子）的3D表示，并可以同时显示分子，这是之前工具所不具备的。这种“所见即所得”的方式从根本上解决了命令行操作“盲人摸象”的痛点。此外，它是一个活跃维护的开源项目，避免了旧工具有的技术栈过时和无法访问的问题。 Q2: 安装和使用PACKMOL-GUI对用户的技术背景有什么要求？ A2: 要求非常低。用户需要预先安装好VMD和PACKMOL。PACKMOL-GUI的安装过程非常简单，只需将下载的文件夹放置到VMD的插件目录中，并在VMD的启动文件中添加一行命令即可。整个过程无需编译，并且有详细的README文件指导。熟悉VMD基本操作的用户可以非常快速地上手。 Q3: 既然PACKMOL-GUI如此强大，它是否存在一些潜在的局限性？ A3: 尽管论文没有专门讨论局限性，但可以推断出几点。首先，它的性能和稳定性完全依赖于VMD。如果VMD在处理超大规模体系（例如数百万原子）时变得卡顿，那么GUI的交互体验也会下降。其次，虽然GUI简化了操作，但正确设置物理化学上合理的约束仍然需要用户的专业知识。例如，在囊泡案例中，如何确定内外层脂质的约束半径，仍然需要用户对手头体系的尺寸有清晰的理解。最后，GUI的最终产物是PACKMOL的输入文件，如果PACKMOL本身在处理某些极端复杂的几何约束时收敛困难，GUI也无法解决这个后端计算的根本问题。关键结论与批判性总结核心结论成功开发了一款名为PACKMOL-GUI的VMD插件，它首次为PACKMOL提供了一个集参数设置、分子可视化和几何约束实时可视化于一体的强大、开源图形用户界面。实现了与VMD的无缝集成，创建了一个直观、高效的一体化平台，用户可以通过“所见即所得”的方式交互式地构建复杂的分子系统。显著提升了建模效率，通过内置的分子数据库、自动分子数估算和清晰的模块化界面，将原本繁琐耗时的命令行操作转变为简单的图形化点击和设置。通过两个复杂的案例研究（球形囊泡和MOF吸附体系），证明了PACKMOL-GUI在处理真实科研问题时的高效性和可靠性。批判性总结与展望 PACKMOL-GUI的出现，无疑是计算化学和分子模拟领域一个极其重要且实用的工程实践成果。它精准地解决了PACKMOL这个“叫好不叫座”（功能强大但使用不便）工具的核心痛点，极大地降低了构建复杂分子体系初始构象的门槛。通过将其巧妙地植入VMD这一事实上的行业标准可视化软件中，作者确保了该工具能被最广泛的科研群体快速接受和使用。可以预见，该插件将极大地促进VMD和PACKMOL的用户群体增长，并成为教授分子模拟课程、进行探索性建模的必备工具。潜在的局限性在于，该工具的价值主要体现在“提效”而非“创新”。它没有改变PACKMOL的算法核心，因此无法解决PACKMOL本身可能存在的收敛性或算法上的难题。未来的发展方向可能包括：1）与更多的分子操纵或模拟设置工具（如Molcontroller的更深度集成）联动，实现更复杂的自动化建模流程。2）引入机器学习模型，根据分子类型和约束形状，智能推荐更优的堆叠策略或参数。3）进一步扩充和维护其内置的分子数据库，使其成为一个更加全面的分子建模资源库。

Molecular Dynamics · 2025-08-15

FE-ToolKit: A Versatile Software Suite for Analysis of High-Dimensional Free Energy Surfaces and Alchemical Free Energy Networks

FE-ToolKit：一个用于分析高维自由能表面和炼金术自由能网络的多功能软件套件 📖 本文基本信息摘要自由能模拟在酶设计、药物发现和生物分子工程等多种生物学应用中发挥着关键作用。要表征复杂酶促反应机理背后的高维自由能表面，需要通过伞形采样或弦方法模拟进行广泛的采样。准确地对大型配体库的靶标结合自由能进行排序，则依赖于组织成热力学网络的全面炼金术自由能计算。这些方法的预测准确性取决于强大且可扩展的工具，用于进行全网络数据分析并从异构模拟数据中提取物理性质。在这里，我们介绍了FE-ToolKit，一个多功能的软件套件，用于自动分析自由能表面、最小自由能路径和炼金术自由能网络（热力学图）。引用信息 Giese, T. J., Snyder, R., Piskulich, Z., Barletta, G. P., Zhang, S., McCarthy, E., Ekesan, Ş., & York, D. M. (2025). FE-ToolKit: A Versatile Software Suite for Analysis of High-Dimensional Free Energy Surfaces and Alchemical Free Energy Networks. Journal of Chemical Information and Modeling, 65(17), 5273–5279. https://doi.org/10.1021/acs.jcim.5c00554 引言在现代计算化学与生物物理学领域，自由能计算是理解和预测分子识别、反应机理及构象动力学等核心科学问题的基石。然而，这些计算本身面临着巨大的挑战，主要源于分子构象空间的广阔性以及对稳健统计方法和严格误差分析的内在需求。为应对这些挑战，FE-ToolKit应运而生。它是一个综合性的集成软件包，旨在为两类主要的计算问题——高维自由能面（Free Energy Surface, FES）的表征和炼金术自由能网络的分析——提供一个模块化、面向工作流程的解决方案生态系统。本报告将遵循该工具包的结构，深入剖析其三个核心组成部分：首先是利用ndfes程序进行自由能面的构建与分析；其次是采用edgembar程序执行可扩展的炼金术网络计算；最后是介绍fetkutils中一系列增强计算效率与数据质量的辅助工具。为了给读者提供一个清晰的概览，下表总结了FE-ToolKit生态系统中的关键程序及其核心功能。表1：FE-ToolKit程序生态系统程序/脚本核心功能 ndfes 使用MBAR/vFEP方法，根据伞形采样数据计算N维FES。 ndfes-path-analyzesims.py 为表面加速弦方法（SASM）提取当前迭代的样本并准备ndfes元文件。 ndfes-path 在静态FES上优化最小自由能路径，并为下一次迭代生成新的模拟输入。 edgembar 对单个炼金术变换（“边”）进行MBAR分析，并生成其有效目标函数。 edgembar-WriteGraphHtml.py 执行炼金术自由能的网络范围分析，并生成交互式HTML报告。 fetkutils-tischedule.py 优化炼金术自由能计算中的λ调度，以提高采样效率。 ndfes-AvgFESs.py 平均多个独立的FES，并根据试验间的方差调整不确定性。 ndfes-CombineMetafiles.py 将多个元文件合并为一个，以聚合采样数据。 ndfes-PrintFES.py 将FES检查点文件中的数据打印为人类可读的文本格式。 Figure 1. FE-ToolKit consists of ndfes for calculating N-dimensional free energy surfaces, edgembar for analyzing alchemical free energy networks using the EdgeMBAR method, and FE-ToolKit utilities (fetkutils) for optimizing schedules of alchemical states. 第一部分：使用ndfes构建和分析自由能形貌 ndfes是FE-ToolKit中用于将偏置模拟（biased simulation）数据转化为有意义的多维自由能面（Free Energy Surface, FES），并识别其上最可能转变路径的核心组件。本节将详细阐述其理论基础、核心方法及实现细节。 1.1 伞形采样与集体变量（CVs）的原理伞形采样（Umbrella Sampling）是一种成熟的增强采样技术，常用于计算分子构象变化、化学反应或分子解离/结合等过程的自由能。它通过施加人工偏置势（biasing potential）来克服高自由能垒，从而确保沿特定过程坐标的充分采样。在FE-ToolKit的语境中，这些过程坐标被称为“反应坐标”（Reaction Coordinates）或更广为人知的“集体变量”（Collective Variables, CVs）。CVs是描述所研究过程的一组低维坐标。本文中的示例并未指定具体的分子体系（例如，某个特定的蛋白或反应），而是作为通用教程进行展示。但其中使用的CVs是该领域的典型代表，例如以埃（Å）为单位的原子间距离，或以度（degrees）为单位的角度或二面角。 1.2 从偏置数据到无偏表面：MBAR与vFEP方法核心问题是如何将来自多个独立的、仅探索了CV空间小范围的偏置模拟数据，整合成一个全局的、无偏的自由能面。FE-ToolKit为此提供了两种功能强大且互为补充的先进方法。多态贝内特接受率（MBAR）方法 MBAR 是一种在统计上被证明是最优的数据重权（reweighting）技术。其核心思想是：所有偏置模拟（每个模拟是一个“态”）的采样数据可以被汇集起来，通过一个优化的权重因子，来估计任何一个“态”（包括我们最感兴趣的无偏物理态）的性质。详细的原理见下一篇推送。变分自由能剖面（vFEP）方法 vFEP是一种参数化方法，与MBAR不同，它不直接计算离散点的概率，而是假设整个自由能表面（FES）可以用一个全局的、连续光滑的解析函数 $F_h(\xi;p)$ 来建模。其核心思想是通过光滑函数拟合能量地貌，类似于用一条平滑的数学曲线或曲面来拟合整个能量地形。 vFEP通过最大化观测到所有偏置模拟样本的对数似然函数来找到最优的函数参数 $p$。该方法使用基数B样条作为基函数来构建全局函数 $F_h(\xi;p)$。B样条是一种标准化的、柔性的“曲线积木”，每个基函数在空间的一小块区域内有值，其他地方为零，特别适合描述规则网格上的函数。详细的原理见下一篇推送。 vFEP与MBAR互补。MBAR是非参数化的，忠实于原始数据，但在数据稀疏区域可能结果不连续或噪声多；vFEP是参数化的，假设FES平滑，能提供平滑连续的表面并便于后续分析，但可能引入模型偏见。用户可根据具体问题选择合适工具或联合使用进行交叉验证。 1.3 寻找最优路径：表面加速弦方法（SASM）在获得了FES之后，下一个重要任务是识别连接两个稳定态（如反应物和产物）的最小自由能路径（Minimum Free Energy Path, MFEP）。ndfes-path 程序为此实现了弦方法的一个先进变体——表面加速弦方法（SASM）。 SASM的迭代过程 SASM的迭代流程针对的是一个特定的反应过程或构象转变（例如一个蛋白的开闭运动，或一个配体的解离路径），而不是一次性处理多个不同的配体。其核心思想是，它将路径（“弦”）的表示与用于生成FES的伞形采样解耦。它的可靠性来源于一个“数据驱动的、渐进精化的”迭代过程：初始猜测与采样：基于一个初始猜测的路径进行初步的伞形采样。聚合与分析：使用ndfes-path-analyzesims.py脚本收集当前及所有先前迭代的全部采样数据。随后，运行ndfes程序，基于这些聚合数据计算出当前对全局FES的最佳估计。路径优化：ndfes-path程序读取步骤2中生成的静态FES，并在此固定的表面上优化弦的位置，以找到当前对MFEP的最佳估计。采样精化与迭代：最后，ndfes-path生成新的模拟输入文件。这些文件会在新优化的路径周围放置新的伞形采样窗口以提高路径局部的分辨率，或在路径的末端进行采样以将其扩展到未探索的区域。随后返回步骤2，进行下一轮迭代。这个策略通过利用全部历史数据来不断修正全局FES，确保路径优化总是在最可靠的表面上进行，从而防止路径在FES定义不清的区域中“迷失”，加速收敛至真实的MFEP。 graph LR %% 定义节点和边的样式 classDef startNode fill:#E8F8F5,stroke:#16A085,stroke-width:2px,font-family:SimHei classDef processNode fill:#EAF2F8,stroke:#5499C7,stroke-width:2px,font-family:SimHei classDef loopArrow stroke:#E74C3C,stroke-width:2.5px,stroke-dasharray: 5 5 %% 节点定义 A("1.初始猜测与采样 基于初始路径进行初步伞形采样") B["2.聚合与分析 程序：ndfes-path-analyzesims.py 与 ndfes 聚合所有历史数据并计算全局FES"] C["3.路径优化 程序：ndfes-path 在固定的FES上优化路径"] D["4.采样精化 程序：ndfes-path 生成新的伞形采样窗口"] %% 流程连接 A --> B B --> C C --> D D -- "返回步骤2 进行下一轮迭代" --> B %% 为节点和边应用样式 class A startNode class B,C,D processNode linkStyle 3 stroke:#c0392b,stroke-width:2px 1.4 最终的自由能面：结构与内容 ndfes的最终输出是一个离散化的多维网格，存储在一个信息详尽的XML格式的检查点文件中。这个输出远不止是能量值，网格中的每个“箱”（bin）都包含了用于分析和质量评估的丰富数据。表1：一个ndfes FES箱的数据结构（MBAR模型）数据字段描述与单位重要性 Bin坐标 (<bidx>) 标识箱在多维网格中位置的一组整数索引。定义了FES上的一个特定离散点。自由能 (<val>) 箱中心的自由能值，单位为kcal/mol。计算的主要结果，描述了该状态的相对稳定性。标准误差 (<err>) 自由能值的不确定性，通过自助法（bootstrap）估计，单位为kcal/mol。衡量结果的统计置信度，是误差分析的关键。 Bin布居数 (<size>) 落入该箱的原始样本数量。表明该区域的采样质量；数量过少可能意味着结果不可靠。重权熵 (<re>) 一个介于0和1之间的无量纲数。衡量不同偏置模拟之间的重叠程度，越接近1越好。第二部分：使用edgembar进行网络范围的炼金术计算 FE-ToolKit 的 edgembar 组件为相对自由能计算提供了一个强大且可扩展的解决方案，尤其适用于处理大规模的配体结合或溶剂化能研究。 2.1 炼金术网络范式为了计算相对结合或溶剂化自由能，通常会构建一个热力学循环。这些计算可以被直观地表示为一个图形网络：网络中的节点（nodes）代表不同的分子（如配体），而连接两个节点的边（edges）则代表在这两个分子之间进行的炼金术转换。每条边关联的值是相对自由能差，记为 $\Delta\Delta G$。例如，在计算相对结合自由能时，该值定义为 $\Delta\Delta G_{(ab)} = \Delta G_{(ab),protein} - \Delta G_{(ab),water}$。这个值直接反映了配体 B 相对于配体 A 与靶标蛋白结合的优势或劣势程度。 2.2 EdgeMBAR 方法：一种可扩展的网络分析策略当我们需要比较一系列配体（例如候选药物分子）与同一靶点的结合能力时，通常会构建一个“炼金术自由能网络”。edgembar 是 FE-ToolKit 中为此类任务量身定做的核心组件。它采用了一种创新性的 EdgeMBAR 方法，将复杂的网络分析问题分解为几个清晰、高效的步骤。 graph TD %% 定义节点样式 classDef inputNode fill:#FEF9E7,stroke:#F39C12,stroke-width:2px,font-family:SimHei classDef stepNode fill:#EAF2F8,stroke:#5499C7,stroke-width:2px,font-family:SimHei classDef innovationNode fill:#E8DAEF,stroke:#8E44AD,stroke-width:2px,font-family:SimHei classDef solveNode fill:#D5F5E3,stroke:#229954,stroke-width:2px,font-family:SimHei classDef outputNode fill:#E8F8F5,stroke:#16A085,stroke-width:2px,font-family:SimHei %% 节点定义 A("炼金术自由能网络 包含所有边的原始模拟数据") B["步骤一：隔离与表征 对每一条边独立进行MBAR分析 得到无约束自由能 g(ab)"] C["步骤二：抽象为有效模型 (核心创新) 将每条边的复杂目标函数 近似为简单的二次函数 提取 g(ab) 和置信度 k(ab)"] D["步骤三：线性代数求解 整合所有边的 g 与 k 信息 构建并求解全局线性方程组"] E("最终网络解 c 得到所有配体全局一致的 有约束自由能(CFE)") %% 流程连接 A --> B B --> C C --> D D --> E %% 为节点应用样式 class A inputNode class B stepNode class C innovationNode class D solveNode class E outputNode 步骤一：隔离与表征（Isolation & Characterization）在这一阶段，edgembar 将复杂的网络拆解开，对其中的每一条“边”（edge）进行独立的、高精度的分析。通俗解释：可以把构建整个自由能网络比作绘制一幅完整的国家地图。传统方法可能试图一次性测量和绘制所有省份，计算量巨大且容易出错。EdgeMBAR 则更像“分而治之”：它先向每个省（每一条边）派出一位独立的“本地勘探专家”。这位专家只负责深度勘探自己省内的地形，完全不受邻省情况的干扰。技术实现：对于网络中任意一条代表“配体 a → 配体 b”转换的边 (ab)，程序首先会构建其完整的“边目标函数” $F_{(ab)}(G_{(ab)})$。该函数是这条边所有相关模拟试验（包括不同环境、阶段和重复试验）的 MBAR 目标函数的总和。通过最小化这个局部的目标函数（$G_{(ab)}^{} = \arg\min F_{(ab)}(G_{(ab)})$），可以得到该边在完全独立、不受网络中其他边影响时的无约束相对自由能（unconstrained relative free energy），记为 $g_{(ab)} = \Delta\Delta G_{(ab)}^{}$。这代表了这条边基于其自身模拟数据得出的“最真实”的自由能值。步骤二：抽象为有效模型（Abstraction to an Effective Model）这是 EdgeMBAR 方法的核心创新所在。在进行全局网络分析时，程序并不直接使用那个包含了所有原始数据、形式复杂的 $F_{(ab)}$，而是用一个极其简单的二次函数（即抛物线）来近似模拟其在最小值 $g_{(ab)}$ 附近的行为。通俗解释：那位“本地勘探专家”在完成详细勘探后，并不会把所有密密麻麻的原始测绘数据都上报给总部。他只提交一份极其凝练的报告：“我省的平均海拔是 $g_{(ab)}$，我对这个值的置信度是 $k_{(ab)}$。” 技术实现：这个近似的二次“有效目标函数”形式如下： \[\tilde{F}_{(ab)}(x) \approx F_{(ab)}(G_{(ab)}^{*}) + \frac{k_{(ab)}}{2}(x - g_{(ab)})^{2}\] 这个简单的抛物线精确地抓住了关于这条边计算结果的两个最关键信息：最可能的自由能值 ($g_{(ab)}$)：即抛物线的最低点位置，代表了独立的边分析给出的最佳估计值。结果的置信度或精度 ($k_{(ab)}$)：这是抛物线的“力常数”，决定了曲线的陡峭程度。$k_{(ab)}$ 越大，抛物线越“瘦削”，意味着微小的自由能偏差 $x - g_{(ab)}$ 都会导致目标函数值急剧上升。这表明模拟数据非常肯定地指向 $g_{(ab)}$ 这个值，因此该计算结果的置信度越高、统计误差越小。反之，一个平坦的抛物线（$k_{(ab)}$ 很小）则代表了较大的不确定性。步骤三：可扩展的线性代数求解（Scalable Linear Algebra Solution）通过将网络中的每一条边都抽象成一个简单的二次函数，原本棘手的、需要处理海量原始数据的非线性优化问题，被神奇地转化为了一个可以解析求解的线性代数问题。通俗解释：总部现在收到了来自所有省份的标准化报告（$g_{(ab)}$ 和 $k_{(ab)}$）。为了绘制一张全局一致的国家地图，总部只需执行一个简单的“加权平均”过程：找到一组能最好地同时满足所有本地报告，且优先采纳那些置信度（$k_{(ab)}$ 值）更高的报告的“官方海拔值”（$c_a, c_b, \dots$）。技术实现：整个网络的全局目标函数现在是所有边的有效目标函数之和，这是一个简单的二次型： \[F_{\text{graph}}(\mathbf{c}) = \sum_{(ab)}^{N_{\text{edge}}} \frac{k_{(ab)}}{2} (c_b - c_a - g_{(ab)})^2\] 其中 $\mathbf{c}$ 是一个包含了所有节点（配体）待求的相对自由能 $c_a, c_b, \dots$ 的向量。最小化这个函数等价于求解一个线性方程组，其闭合解形式非常简洁： \[\mathbf{c} = \mathbf{M}^{-1} \cdot \mathbf{X}^T \cdot \mathbf{K} \cdot \mathbf{g}\] 这里的 $\mathbf{g}$ 是所有无约束自由能构成的向量，$\mathbf{K}$ 是所有力常数构成的对角矩阵，$\mathbf{X}$ 和 $\mathbf{M}$ 是描述网络拓扑结构（即节点如何被边连接）的矩阵。这种方法的优势是巨大的：计算效率极高：求解线性方程组远比处理海量原始数据和最小化非线性函数要快得多。出色的可扩展性：如果网络中增加了一条新的边，我们只需对这条新边执行一次步骤一和步骤二，然后几乎可以瞬时解出新的全局网络结果。而传统方法可能需要从头重新分析整个网络，成本高昂。步骤四：得到网络解以后能做什么？——从无约束到约束分析求解出向量 $\mathbf{c}$（即所有配体的相对自由能）后，我们可以进行一系列深刻的分析，这正是edgembar的核心价值所在。计算有约束自由能（Constrained Free Energy, CFE）求解该线性方程组的主要目的之一就是计算有约束自由能。定义：向量 $\mathbf{c}$ 中的解，代表了在满足全局热力学循环闭合条件下，对所有配体相对自由能的最佳估计。网络中任意两点（例如配体 a 和配体 b）的自由能差 $c_b - c_a$，就是该路径的有约束自由能（CFE）。与无约束自由能（UFE）的对比：与之对应的是我们在步骤一中得到的无约束自由能（UFE），即 $g_{(ab)}$。UFE 是单条边“认为”自己应该有的值，而 CFE 是它在整个关系网中为了“合群”（满足热力学一致性）而必须调整到的值。诊断价值：比较 CFE 和 UFE 的差异，即 Shift（$ CFE - UFE $），是一个极其重要的诊断指标。一个很大的 Shift 值意味着这条边的独立计算结果与网络中的其他邻居存在严重冲突，表明这条边的模拟数据可能存在问题，需要仔细检查。整合实验数据进行进一步约束 edgembar 的强大之处还在于，它允许将外部的高精度数据（如已知的实验测量值）作为额外的、更强的约束条件整合到网络分析中。实现机制：该功能通过拉格朗日乘子法（Lagrange’s method of undetermined multipliers）实现。它在最小化全局目标函数 $F_{\text{graph}}(\mathbf{c})$ 的同时，额外施加了一组线性约束，强制要求网络中某些边的 CFE 值必须等于给定的实验值。实际操作：用户可以通过在 edgembar-WriteGraphHtml.py 脚本中使用 --expt FILENAME 和 --constrain LIGA~LIGB 等命令行选项来轻松实现这一功能。意义：这使得我们可以利用已知的、可靠的实验数据来“锚定”整个自由能网络，从而提高对未知配体自由能的预测精度。深入的诊断与质量评估最终的“graph report”（HTML 格式的图报告）提供了丰富的诊断指标，帮助用户全面评估网络质量。表 2：网络分析中的关键诊断指标指标全称描述与意义 UFE / dUFE Unconstrained Free Energy 边的无约束自由能及其标准误，来自独立的边分析。 CFE / dCFE Constrained Free Energy 边的有约束自由能及其标准误，来自网络全局分析的结果。 Force Constant ($k_{(ab)}$) 有效目标函数中的二次项系数，反映了自由能曲线的陡峭程度。 “力常数”越小，表示独立计算该边自由能的不确定性越大，其结果在网络整合中的权重也相应较低。 Shift Shift 网络自洽后的边自由能与独立计算的边自由能之差的绝对值：$ \Delta\Delta G_{\text{network}} - \Delta\Delta G_{\text{isolated}} $。该值较大时，表明网络整合显著改变了该边的自由能估计，可能暗示网络中存在不一致性或该边的初始计算存在偏差。 CC Cycle Closure error 任何一个闭合环路的 UFE 之和的绝对值，直接衡量网络的不一致程度。 Average Cycle Closure (ACC) 遍历某条边的所有独立闭合路径的循环闭合误差的平均值。 ACC 值较大同样标志着该边是网络不一致性的主要来源，需要仔细检查与之相关的模拟数据。 LMI Lagrange Multiplier Index 一个无量纲数，衡量一条边对整个网络施加的“应力”或“张力” 。值越大，表明该边与网络其余部分的矛盾越大。 OFC2 Objective Force Constant 目标函数力常数的两倍 (2k(ab))，衡量 UFE 计算结果的置信度。 2.3 实用的网络分析与诊断 FE-ToolKit 的设计理念是赋能专家用户，因此 edgembar 及其配套脚本不仅提供最终的自由能数值，还输出了大量的诊断数据，以评估结果的可靠性和整个网络的一致性。edgembar-WriteGraphHtml.py 脚本生成的交互式 HTML 报告是一个强大的可视化工具，用户可以用它来探索网络图、节点、边和循环的详细属性。为了有效利用这些诊断信息，理解关键指标的含义至关重要。补充细节：edgembar 的输入与输出输入要求： edgembar 的输入是一个 XML 文件，该文件组织模拟数据到环境、阶段、试验和状态的层次结构中。每个试验的数据是一组名为 “efep_tlam_elam.dat” 的文件集合，其中 tlam 是采样状态，elam 是文件中制表的势能状态。文件的第一列是模拟时间（皮秒），第二列是势能（kcal/mol）。如果需要，还可以包含额外列用于不同环境和目标势能。输出与报告： edgembar 的输出被组织成数据结构并写入 Python 文件，可直接导入其他脚本进行分析。执行 Python 输出会总结结果到一个 HTML 格式的 “边报告” 中。 edgembar-WriteGraphHtml.py 脚本读取多个 edgembar 输出，计算配体自由能，并总结结果到一个 HTML 格式的 “图报告” 中，比较孤立边自由能与配体自由能差异，并包含闭合路径及其闭合误差的表格。实际应用案例假设我们正在进行一项大规模的配体结合自由能计算，以筛选潜在的药物分子。我们构建了一个包含 100 个配体的网络，每个配体与相邻配体之间都有边连接，形成一个复杂的热力学网络。使用 edgembar，我们可以：对每条边进行独立分析，计算其无约束相对自由能 $g_{(ab)}$ 和力常数 $k_{(ab)}$。将每条边的结果抽象为二次有效目标函数，构建整个网络的全局目标函数。求解线性方程组，得到所有配体的相对自由能。利用 edgembar-WriteGraphHtml.py 生成的 HTML 报告，检查 Shift、LMI 和 ACC 等诊断指标，识别网络中的潜在问题边。针对问题边进行进一步的模拟或调整计算参数，优化网络一致性。通过这种系统性的分析和诊断流程，edgembar 不仅提供了准确的相对自由能计算结果，还帮助研究人员深入理解网络中各边和节点的相互作用，为药物设计和分子模拟提供了宝贵的指导。第三部分：辅助工具与实用工作流程 (fetkutils) 如果说 ndfes 和 edgembar 是执行核心科学分析的“主力部队”，那么 fetkutils 工具集就是保障整个研究工作流程顺畅、高效、可靠的“精英后勤与工程团队”。它解决了两个在自由能计算中普遍存在、至关重要的实践问题：如何用最少的资源达到最高的计算效率，以及如何确保用于分析的数据是稳定可靠的。优化模拟效率：“智能的领航员” (fetkutils-tischedule.py) 核心思想：与其亡羊补牢，不如未雨绸缪。在进行昂贵的炼金术自由能计算时，一个常见的效率瓶颈是不同炼金术状态（λ态）之间的转换不顺畅。可以把这个过程想象成一场长距离接力赛，如果其中某两个赛段的交接棒（状态交换）非常困难，那么整个团队的速度都会被拖慢。天真地将“接力点”（λ值）均匀分布，往往不是最高效的策略。 fetkutils-tischedule.py 工具提供了一种主动优化的智能策略。它就像一位经验丰富的教练，在正式比赛前，先让队员们进行一次简短的“测试跑”（即“预烧”模拟），以识别出哪些交接棒环节是薄弱点。然后，它利用这些测试数据，为正式比赛量身定做一套最优的接力方案（即优化的λ调度表）——在困难的交接区段，让接力点靠得更近，在轻松的区段则拉得更远。这个“先侦察，再规划”的策略，能够确保最终进行的、计算成本高昂的生产性模拟从一开始就在最优化的路径上运行，从而显著节省宝贵的计算资源和研究时间。确保数据质量：“严谨的质检员”（自动化平衡检测）核心思想：用客观的算法取代主观的人眼判断。分子模拟的轨迹数据，就像刚从工厂生产出来的产品，必须经过严格的质量检验才能使用。每条轨迹的开头部分都是系统从初始状态走向平衡的“预热”或“适应”阶段，这部分数据是不稳定、不可靠的，必须被准确地切除。在面对成百上千条模拟轨迹时，手动检查并决定切割点不仅繁琐，而且极易引入研究者的主观偏见。 FE-ToolKit 内置的自动化平衡检测算法就是一位不知疲倦且铁面无私的“质检员”。它会自动审查每一条轨迹的关键数据流（如能量波动），并运用一套严格的统计检验流程来做出判断。它会反复“考察”数据，直到找到一个明确的、不再有系统性漂移或剧烈波动的稳定“生产区域”。这个自动化流程提供了一种可重复的、客观的方法来筛选数据，从源头上保证了只有高质量、已平衡的样本才能被用于最终的自由能分析，这对于确保科学结论的可靠性至关重要。其他工具功能概览 ndfes-AvgFESs.py：用于合并来自独立重复试验的结果，并正确地传递误差，这对于评估结果的稳健性至关重要。 ndfes-CombineMetafiles.py：一个实用的工具，用于聚合来自多个模拟集的数据，简化了对大规模伞形采样数据的管理。 ndfes-PrintFES.py：用于从二进制的检查点文件中提取数据，并将其转换为人类可读的文本格式，方便后续处理或绘图。 ndfes-genbias：一个专门用于处理非谐波偏置势的程序。这体现了工具包的灵活性，但文档也明确指出，使用该程序会在效率和易用性上有所取舍。它们共同构成了面向工作流程的完整工具链。结论 FE-ToolKit不仅仅是一个程序的集合，它体现了对现代自由能计算的一种连贯而强大的构想。通过对其核心组件和设计理念的深入剖析，可以总结出几个贯穿始终的主题：可扩展的严谨性：无论是通过SASM中的解耦策略，还是EdgeMBAR中革命性的抽象方法，该工具包始终在追求统计上最优的严谨性的同时，确保方法能够扩展到日益复杂的系统中。赋能专家用户：从提供MBAR和vFEP两种FES构建方法，到输出详尽的网络诊断指标，FE-ToolKit的设计处处体现了对专业用户的尊重，为他们提供了深入分析和验证计算结果所需的全部工具。模块化的工作流程设计：工具包由一系列目标明确、可协同工作的脚本和程序组成，形成了一个从实验设计（如λ调度优化）、数据生成、核心分析到最终结果可视化的完整生态系统。抽象的力量：EdgeMBAR方法是这一点的最佳体现。通过将复杂的边目标函数抽象为一个简单的二次模型，它成功地将一个难以处理的全局优化问题转化为一个易于求解的线性问题，这正是计算科学中优雅解决方案的典范。综上所述，FE-ToolKit为计算科学家提供了一个从头至尾的解决方案，引导研究人员高效、自信地应对从基础反应机理到大规模药物设计等领域中极具挑战性的自由能计算问题。局限性与未来展望尽管FE-ToolKit功能强大，但根据原文的描述，其在当前版本中仍存在一些局限性，并指明了未来的发展方向：特定组件的功能限制：工具包中提供了一个用于处理通用偏置势的程序 ndfes-genbias，但作者明确建议除非绝对必要，否则不推荐使用。主要原因是 ndfes-genbias 需要写入非常大的输入文件，对内存工作站的要求很高。此外，该程序尚不能执行vFEP方法，并且在聚合来自多个重复试验的数据时，由于“偏置索引”可能会失效，操作起来非常谨慎和困难。性能与实现：核心的网络分析程序 edgembar 是一个用C++编写的、支持OpenMP并行的实现，但原文明确指出它目前缺乏GPU加速功能。在当前大规模计算日益依赖GPU的背景下，这可能在处理超大规模网络时成为一个潜在的性能瓶颈。当前版本的功能待完善之处：对于在不同温度下进行的模拟，并试图在某个特定温度下分析其自由能面的功能，原文提到这部分功能尚未经过广泛测试，且初步测试表明结果可能会受到显著的数值噪声影响。在能量单位方面，当前版本的图报告和边报告中的能量单位是固定的（kcal/mol）。原文提到未来的版本将允许用户选择输出的能量单位，暗示了当前版本在这方面缺乏灵活性。持续发展的需求：作者在结尾处指出，FE-ToolKit将继续被开发和维护，以支持新兴的集成自由能方法。这表明该工具包虽然在处理当前主流方法上非常成熟，但仍需不断迭代，以跟上计算化学领域快速涌现的新技术和新方法。

Molecular Dynamics · 2025-06-23

FE-ToolKit Methodology Deep Analysis: Mathematical Derivations and Physical Interpretations

Molecular Dynamics · 2025-06-23

GROMACS Defaults in Topology Files: Understanding comb-rule and nonbond_params

GROMACS 中 comb-rule 与 [nonbond_params] 参数解析本文档旨在详细解释 GROMACS 拓扑文件中 [defaults] 指令下的 comb-rule（组合规则）以及 [atomtypes] 和 [nonbond_params] 部分中非键参数（特别是 Lennard-Jones 参数）的含义和解释方式。一、[defaults] 指令详解在 GROMACS 的拓扑文件（通常是 .top 文件或力场主 .itp 文件）中，[defaults] 指令用于设定非键相互作用的全局默认行为。示例 [ defaults ] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 2 no 1.0 1.0 参数解释 nbfunc (Non-bonded function type)：定义非键势函数类型。 1：Lennard-Jones 势。这是绝大多数经典力场（如 AMBER, CHARMM, OPLS, Martini）使用的形式。 2：Buckingham 势。注意：根据 GROMACS 文档和社区讨论，Buckingham 势 (nbfunc = 2) 自 GROMACS 2019 版本后可能已被弃用或不再完全支持。参考链接：https://gromacs.bioexcel.eu/t/how-use-desired-mixing-rule-in-gromacs/10409/3 comb-rule (Combination rule)：定义当 [nonbond_params] 部分没有显式给出不同原子类型 i 和 j 之间的非键参数时，如何从各自的原子类型参数（[atomtypes] 部分的参数）计算出交叉项参数。 gen-pairs (Generate 1-4 pairs)：决定是否自动生成1-4相互作用对（即通过3个键连接的原子对）。 yes：根据成键信息自动生成，并通常与 fudgeLJ 和 fudgeQQ 联用。 no：不自动生成，1-4相互作用需要在 [pairs] 或 [pairtypes] 部分显式定义，或者由力场设计本身通过其他方式处理（如Martini）。 fudgeLJ：如果 gen-pairs = yes，此参数定义了1-4相互作用中 Lennard-Jones 部分的缩放因子。 fudgeQQ：如果 gen-pairs = yes，此参数定义了1-4相互作用中静电部分的缩放因子。 GROMACS comb-rule：对 [atomtypes] 参数的解释及交叉项的计算 comb-rule 的设置直接影响 GROMACS 如何解释 [atomtypes] 部分中的 V 和 W 列参数，以及在 [nonbond_params] 中没有显式定义一对原子类型间的非键参数时，如何计算这些交叉项参数。 https://manual.gromacs.org/current/reference-manual/topologies/parameter-files.html#non-bonded-parameters 1. [atomtypes] 中 V 和 W 参数的解释根据 GROMACS 手册: 如果 comb-rule = 1: $V_{ii}$ 代表 $C_{6,ii} = 4 \epsilon_{ii} \sigma_{ii}^6$ (单位：kJ mol⁻¹ nm⁶) $W_{ii}$ 代表 $C_{12,ii} = 4 \epsilon_{ii} \sigma_{ii}^{12}$ (单位：kJ mol⁻¹ nm¹²) 此时 Lennard-Jones 势能通常写作： \[V_{LJ}(r) = \frac{C_{12,ij}}{r^{12}} - \frac{C_{6,ij}}{r^6}\] 如果 comb-rule = 2 或 3: $V_{ii}$ 直接代表 $\sigma_{ii}$ (单位：nm) $W_{ii}$ 直接代表 $\epsilon_{ii}$ (单位：kJ mol⁻¹) 此时 Lennard-Jones 势能通常写作： \[V_{LJ}(r) = 4 \epsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r}\right)^{12} - \left(\frac{\sigma_{ij}}{r}\right)^6\right]\] 2. 交叉项参数的计算 (如果未在 [nonbond_params] 中显式定义) 对于 comb-rule = 1 和 3: GROMACS 使用几何平均来组合 $C_6$ 和 $C_{12}$ 参数： \[C_{6,ij} = \sqrt{C_{6,ii} \times C_{6,jj}}\] \[C_{12,ij} = \sqrt{C_{12,ii} \times C_{12,jj}}\] 注意：如果 comb-rule = 3，[atomtypes] 中的 $V_{ii}$ 和 $W_{ii}$ 被解释为 $\sigma_{ii}$ 和 $\epsilon_{ii}$。GROMACS 内部会先将它们转换为 $C_{6,ii}$ 和 $C_{12,ii}$，然后再应用上述几何平均规则。对于 comb-rule = 2 (Lorentz-Berthelot 规则): GROMACS 使用算术平均组合 $\sigma$ 参数，使用几何平均组合 $\epsilon$ 参数： \[\sigma_{ij} = \frac{\sigma_{ii} + \sigma_{jj}}{2}\] \[\epsilon_{ij} = \sqrt{\epsilon_{ii} \times \epsilon_{jj}}\] 关于常见力场的组合规则说明注意：常见力场（CHARMM、AMBER、OPLS等）与 GROMACS 中 comb-rule 参数的对应关系在文献中并不十分明晰，以下信息基于有限的资料整理推测：力场 σ 组合规则 ε 组合规则可能的 GROMACS 设置备注 CHARMM 算术平均几何平均 comb-rule = 2 如果 [atomtypes] 中提供的是 $\sigma_{ii}$ 和 $\epsilon_{ii}$ AMBER 算术平均几何平均 comb-rule = 2 明确使用 Lorentz-Berthelot 规则 OPLS 几何平均几何平均 comb-rule = 3 通常在 [nonbond_params] 中显式定义所有交叉项算术平均是Lorentz提出的，几何平均是Berthelot提出的也就是说，comb-rule = 1当然是万能的，但全原子一般是给出 $\sigma$ 和 $\epsilon$，其中comb-rule = 2 即Lorentz-Berthelot 规则，comb-rule = 3 即均为几何平均。 CHARMM：使用 Lorentz-Berthelot 规则。对 $\sigma$ (或NAMD里面，等效的 $R_{min}$) 使用算术平均，对 $\epsilon$ 使用几何平均。 $R_{min,ij} = \frac{R_{min,ii} + R_{min,jj}}{2} \text{(等效于 $\sigma$ 的算术平均)}$ \[\epsilon_{ij} = \sqrt{\epsilon_{ii} \times \epsilon_{jj}}\] 参考：NAMD Mailing List - https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2009-2010/3885.html AMBER：明确使用 Lorentz-Berthelot 规则。根据 AMBER 手册节选：”For Amber force fields, cross terms involving different atom types i and j are evaluated according to the Lorentz/Berthelot mixing rules…“，可以自行查找 OPLS：OPLS 力场通常对 Lennard-Jones 参数 $\sigma$ 和 $\epsilon$ 都使用几何平均。 OPLS 力场的 GROMACS 实现通常没有 [nonbond_params] 。 NAMD参考文末二、[atomtypes] 和 [nonbond_params] 中的参数解释 GROMACS 通过 [atomtypes] 和 [nonbond_params] (或 [pairtypes]) 这两个主要部分来定义非键相互作用参数。参考：GROMACS Manual - Non-bonded parameters - https://manual.gromacs.org/current/reference-manual/topologies/parameter-files.html#non-bonded-parameters [atomtypes] 部分此部分定义了每种原子类型自身 (ii) 的基本非键参数。这些参数的解释（是 $\sigma$, $\epsilon$ 还是 $C_6$, $C_{12}$）取决于 [defaults] 中设置的 comb-rule。示例 (OPLS-AA 风格，通常 comb-rule = 1，意味着 V, W 是 $C_6$, $C_{12}$) [ atomtypes ] ;name at.num mass charge ptype V(c6) W(c12) ; V 和 W 的含义取决于 comb-rule O 8 15.99940 0.000 A 0.22617E-02 0.74158E-06 ; V(c6) = C6_ii, W(c12) = C12_ii OM 8 15.99940 0.000 A 0.22617E-02 0.74158E-06 ... [nonbond_params] 部分此部分用于显式定义特定原子类型对 i 和 j 之间的非键相互作用参数。这里定义的参数将覆盖任何通过组合规则计算得到的参数。示例1 ( comb-rule = 1 配合，参数为直接的 $C_{6,ij}$ 和 $C_{12,ij}$) [ nonbond_params ] ; i j func V(c6) W(c12) ; 列标题指明了是 C6 和 C12 O O 1 0.22617E-02 0.74158E-06 ; O-O 相互作用的 C6_ij 和 C12_ij O OA 1 0.22617E-02 0.13807E-05 ; O-OA 相互作用的 C6_ij 和 C12_ij ... V(c6)：该原子类型对的 Lennard-Jones $C_{6,ij}$ 参数 (单位：kJ mol⁻¹ nm⁶)。 W(c12)：该原子类型对的 Lennard-Jones $C_{12,ij}$ 参数 (单位：kJ mol⁻¹ nm¹²)。示例2 (Martini 风格，参数为直接的 $\sigma_{ij}$ 和 $\epsilon_{ij}$) [ nonbond_params ] ; i j func sigma epsilon ; 列标题通常会指明是 sigma 和 epsilon P6 P6 1 0.470 4.990 ; P6-P6 相互作用的 sigma_ij 和 epsilon_ij P6 P5 1 0.470 4.730 ; P6-P5 相互作用的 sigma_ij 和 epsilon_ij ... i, j：相互作用的原子类型。 func：函数类型，1 表示 Lennard-Jones 12-6 势。 sigma：该原子类型对的 Lennard-Jones $\sigma_{ij}$ 参数 (单位：nm)。 epsilon：该原子类型对的 Lennard-Jones $\epsilon_{ij}$ 参数 (单位：kJ/mol)。关键点：[nonbond_params] 中参数的含义（是 $\sigma$/$\epsilon$ 还是 $C_6$/$C_{12}$）直接由该力场文件在该部分的列定义（通常通过注释中的列标题）决定。func=1 只是表示它是一个12-6型的Lennard-Jones势，但参数的表达形式可以有两种。三、Martini 力场的特殊性对于 Martini 力场 (例如 martini_v3.0.0.itp)：参考文献：PCT Souza, et al., Nat. Methods, 2021. DOI：10.1038/s41592-021-01098-3 （看SI的表） [defaults] 指令 Martini 3 的主 .itp 文件通常包含： [ defaults ] ; nbfunc comb-rule 1 2 ; (通常 gen-pairs no, fudgeLJ/QQ 不适用或设为1.0) 这里的 comb-rule = 2 设定了默认的参数类型。 [atomtypes] 部分（真实示例）在 Martini 3 中，[atomtypes] 部分的 $\sigma$ 和 $\epsilon$ 值都设为 0.0，因为 Martini 的核心在于珠子类型之间的相互作用矩阵： [ atomtypes ] ; name mass charge ptype sigma epsilon P6 72.0 0.000 A 0.0 0.0 P5 72.0 0.000 A 0.0 0.0 ... 这里的 sigma 和 epsilon 都是 0.0，表明它们仅是占位符。 [nonbond_params] 部分（真实示例）这是 Martini 力场定义非键相互作用的关键。Martini 不依赖 GROMACS 的组合规则来生成不同珠子类型之间的相互作用参数。相反，它在 [nonbond_params] 部分显式地定义每一对珠子类型之间的 $\sigma_{ij}$ 和 $\epsilon_{ij}$： [ nonbond_params ] P6 P6 1 4.700000e-01 4.990000e+00 P6 P5 1 4.700000e-01 4.730000e+00 ... 注意这里：没有列标题注释，但根据 Martini 文档，这些参数是 $\sigma_{ij}$ (第4列) 和 $\epsilon_{ij}$ (第5列) 所有珠子对的相互作用都被显式定义因此，当 grompp 处理 Martini 拓扑时，它会优先使用 [nonbond_params] 中为特定珠子对定义的 $\sigma_{ij}$ 和 $\epsilon_{ij}$。只有当某一对珠子类型的相互作用没有在 [nonbond_params] 中显式定义时，才会退回到使用 [defaults] 中指定的 comb-rule 和 [atomtypes] 中的参数来尝试计算（但由于 [atomtypes] 中的值都是 0.0，实际上不会产生有意义的相互作用）。详见上一篇：总结对于标准的 Martini 3 力场文件： [atomtypes] 中的 $\sigma$/$\epsilon$ 都是 0.0：它们是占位符，不用于计算。核心的异类珠子对相互作用参数来自 [nonbond_params]：这是Martini设计的核心。 [nonbond_params] 中提供的是针对特定珠子对 ij 的 $\sigma_{ij}$ 和 $\epsilon_{ij}$：这些不是 $C_{6,ij}$ 和 $C_{12,ij}$。 [defaults] 中的 comb-rule = 2 在 Martini 中更多的是一个形式上的设定：因为所有相关的珠子对相互作用参数都是在 [nonbond_params] 中显式提供的。四、总结：如何判断参数类型判断 .itp 文件中非键参数是 ($\sigma$, $\epsilon$) 还是 ($C_6$, $C_{12}$) 的关键步骤： 1. 查看 [defaults] 指令中的 comb-rule 如果 comb-rule = 1，那么 [atomtypes] 中的 V 和 W 列倾向于是 $C_{6,ii}$ 和 $C_{12,ii}$。如果 comb-rule = 2 或 3，那么 [atomtypes] 中的 V 和 W 列倾向于是 $\sigma_{ii}$ 和 $\epsilon_{ii}$。 2. 仔细阅读 [atomtypes] 和 [nonbond_params] 部分的列标题注释如果列标题明确写着 sigma 和 epsilon，那么这些值就是 $\sigma$ 和 $\epsilon$。如果列标题明确写着 V(c6) 和 W(c12)，那么这些值就是 $C_6$ 和 $C_{12}$。假定开发者不至于搞错，这是最直接的判断依据。 3. 查阅相应力场的原始文献和手册这是最权威的判断依据。力场开发者会明确说明其参数的定义和使用方式。实用建议对于您的脚本而言，如果它需要同时处理可能来自不同力场的 .itp 文件，建议：通过一个参数来指定当前处理的ITP文件中的非键参数是哪种类型或者通过智能解析列标题来判断对于 Martini 这样的特殊情况（[atomtypes] 中都是 0.0），直接使用 [nonbond_params] 中的参数其他参考资料 GROMACS Manual - MDP Options for LJ-PME combination rule: https://manual.gromacs.org/current/user-guide/mdp-options.html#mdp-lj-pme-comb-rule NAMD 的相关设置 NAMD Mailing List：https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2009-2010/3885.html “Yes, as is standard for the CHARMM force field NAMD uses arithmetic mean for sigma, geometric mean for epsilon by default. You can change this by adding ‘vdwGeometricSigma yes’ in the config file to support, e.g., OPLS.” NAMD User Guide：https://www.ks.uiuc.edu/Research/namd/3.0.1/ug/node25.html#7012 “vdwGeometricSigma：Use geometric mean, as required by OPLS, rather than traditional arithmetic mean when combining Lennard-Jones sigma parameters for different atom types.”

Molecular Dynamics · 2025-05-31

Martini 3 Bead Types and Naming Conventions: A Comprehensive Guide

主要参考资料： https://doi.org/10.1038/s41592-021-01098-3 Supporting information for: Martini 3: A General Purpose Force Field for Coarse-Grained Molecular Dynamics https://github.com/ricalessandri/Martini3-small-molecules/blob/main/tutorials/building_block_table.pdf https://advanced.onlinelibrary.wiley.com/doi/full/10.1002/adts.202100391 https://cgmartini.nl/docs/tutorials/Martini3/Small_Molecule_Parametrization/ 1. 引言 (Introduction) Martini 力场是一种广泛应用于生物分子模拟的粗粒化 (Coarse-Grained, CG) 模型 (1)。近年来，经过大幅改进和重新参数化的 Martini 3 版本正式发布 (1)。Martini 3 旨在提供一个通用性更强的 CG 力场，不仅适用于脂质、蛋白质、核酸和糖类等生物大分子体系 (4)，也拓展到了对多种小分子、碳纳米材料以及聚合物的研究 (7)。相较于早期版本，Martini 3 的核心改进包括更优化的非键相互作用平衡、引入了新的珠子 (bead) 类型（包括不同尺寸和化学特性的珠子）、并增强了对特定相互作用（如氢键和电子极化效应）的描述能力 (1)。这些改进使得 Martini 3 能够更准确地预测分子的堆积模式和相互作用，从而在更广泛的应用领域中提供可靠的模拟结果 (1)。Martini 模型通常采用“四对一”的映射方案，即平均四个重原子及其相连的氢原子被粗粒化为一个相互作用中心（珠子），但对于环状结构等特殊化学基团，也支持更高分辨率的映射 (2)。本报告旨在详细阐述 Martini 3 力场中珠子的类型、命名方式的传统和原则，并深入探讨其参数化策略和分子映射方法。最后，将通过一个具体的聚合物——聚[2-(N-氧化-N,N-二乙基氨基)甲基丙烯酸乙酯] (poly[2-(N-oxide-N,N-diethylamino)ethyl methacrylate])——的映射实例，展示如何将这些理论知识应用于实践。 2. Martini 3 核心珠子 (Bead) 类型与命名传统 (Martini 3 Core Bead Types and Naming Conventions) Martini 3 模型的基石在于其多样化的珠子类型，这些珠子代表了不同化学性质的分子片段。理解这些珠子的分类、尺寸和命名规则对于正确构建和解读 CG 模型至关重要。 2.1 主要珠子类型 (Main Bead Types) 与早期版本类似，Martini 3 保留了基于极性的四种主要珠子类型 (8)： P (Polar): 极性珠子，代表强极性基团。 N (Non-polar/Intermediate polarity): 中等极性或非极性珠子，代表具有一定极性或非极性的基团。 C (Apolar/Carbon-like): 非极性珠子，通常代表疏水性的烷烃链等。 Q (Charged): 带电荷珠子，代表离子化的基团。在 Martini 3 中，这些主要类型得到了扩展和细化，引入了新的专用珠子类型 (8)： W (Water): 特定的水珠子，与 Martini 2 中的极性 P4 珠子不同，W 珠子经过独立优化，避免了旧模型中水在室温下结冰等问题。 D (Divalent ions): 二价离子珠子。 X (Halo-compounds): 用于描述含卤素原子的基团。这些主要类型（P, N, C, Q, X）进一步划分为多个亚型，通过数字后缀（通常为1-6）表示其相对极性程度或相互作用强度，数字越大通常表示极性越强或相互作用越强 (10)。Martini 3 将可能的相互作用能级从 Martini 2 的10个扩展到了22个，从而可以更精细地描述不同化学基团间的相互作用差异 (8)。此外，还引入了字母后缀来表征特定的化学特性，如氢键给体/受体能力和电子极化效应 (8)。 2.2 珠子尺寸 (Bead Sizes) Martini 3 引入了三种不同尺寸的珠子，以适应不同分辨率的粗粒化需求，这对于精确描述分子形状和堆积至关重要 (8)： Regular (R): 常规尺寸珠子，其 Lennard-Jones (LJ) 参数中的$\sigma$值约为 0.47 nm。设计用于标准的“4对1”映射，即代表4个重原子及其相连的氢原子。 Small (S): 小尺寸珠子，$\sigma$值约为 0.41 nm。设计用于“3对1”的映射，即代表3个重原子。 Tiny (T): 微小尺寸珠子，$\sigma$值约为 0.34 nm。设计用于“2对1”的映射，即代表2个重原子。这三种尺寸的珠子之间的交叉相互作用 (R-S, R-T, S-T) 都经过了专门的参数化，以确保整个力场的平衡性 (8)。小尺寸 (S) 和微小尺寸 (T) 珠子特别适用于描述环状结构（如芳香环和脂肪环）以及其他需要更高分辨率的线性或支链化学基团 (4)。对于完全支化的片段（如季碳原子或叔胺基团），如果片段包含四个非氢原子，通常会使用较小的珠子（如 S 型珠子），因为中心原子的环境暴露程度降低，其对整体相互作用的影响也相应减小 (8)。 2.3 命名约定 (Naming Conventions) Martini 3 珠子的命名遵循一套系统的规则，通常结合了其尺寸、基本化学类型、极性水平以及特殊功能：尺寸前缀: R: 代表常规尺寸 (Regular)，但在很多情况下，如果珠子名称没有明确的尺寸前缀，则默认为常规尺寸。 S: 代表小尺寸 (Small)。 T: 代表微小尺寸 (Tiny)。基本类型字母: P, N, C, Q, X, W, D，如上所述。极性/相互作用能级: 通常是一个数字（1到6，对于P, N, C, Q, X 类型），表示相对极性或相互作用强度。例如，P1 表示低极性极性珠子，P6 表示高极性极性珠子 (10)。氢键后缀: 用于描述珠子的氢键能力 (10)。 d (donor): 表示氢键给体。 a (acceptor): 表示氢键受体。 da: 表示同时具有氢键给体和受体能力。 0 (zero): 对于Q类型珠子 (如 Q0)，表示不具有特定的氢键给体或受体能力。对于P和N类型珠子，若无 ‘d’ 或 ‘a’ 后缀，通常意味着其氢键能力不是其主要特征，或作为一般极性珠子处理。电子极化效应后缀: 主要用于 C 和 X 类型珠子，以模拟邻近化学基团的诱导/共轭效应对分子片段相互作用的影响，并能捕捉优先取向和卤键能力 (8)。 e (electron-donor/enriched): 表示富电子特性。 v (electron-acceptor/vacancy): 表示缺电子特性。例如，萘中心的珠子类型为 TC5e，表示这是一个富电子的微小尺寸非极性珠子。其他特殊后缀: h: 在某些特定珠子类型中使用，例如在脂质尾链中，C5h 和 C4h 用来区分包含不同数量双键的片段 (12)。 r: 在某些溶剂模型中出现，如甲醇 (MEOH) 用 SP2r 表示 (13)。一个典型的 Martini 3 珠子名称组合了这些元素，例如：SP1d (小尺寸、极性类型、1级极性、氢键给体)，TC5e (微小尺寸、非极性类型、5级相互作用、富电子)。 2.4 Martini 3 珠子类型汇总表 (Comprehensive Table of Martini 3 Bead Types) 为了更清晰地展示 Martini 3 中常用珠子的特性，下表总结了部分代表性珠子类型及其关键属性和通常代表的化学片段。此表并非详尽无遗，更完整的列表和特定分子的参数化可以在 Martini 官方网站和相关出版物中找到 (10)。构建新分子模型时，应参考最新的官方 martini_v3.0.0.itp 文件和相关文献中的构建模块表 (8)。珠子名称 (Bead Name) 主要类型 (Main Type) 尺寸 (Size) 极性水平 (Polarity Level) 氢键 (H-bond) 其他后缀 (Other Suffix) 典型化学基团/片段 (Typical Chemical Group/Fragment) W W R N/A N/A 水 (代表4个水分子) TW W T N/A N/A 微小水 (代表2个水分子)，用于受限空间 C1 C R 1 None 饱和烷烃片段 (-CH2-CH2-CH2-CH2-) SC3 C S 3 None 脂肪环片段 (如环己烷中的 -CH2-CH2-CH2- 单元)，支链烷烃 TC5 C T 5 None 芳香环中的 -CH=CH- 片段 (如苯)，共轭体系 TC5e C T 5 None e 富电子芳香片段 (如萘的中心) P1 P R 1 Donor/Acceptor 弱极性基团，如醚 (-O-) SP2d P S 2 Donor 中等极性氢键给体，如伯醇 (-CH2OH 中的 OH 部分，若映射为S珠) TP4a P T 4 Acceptor 强极性氢键受体，如羰基 (C=O，若映射为T珠) N0 N R 0 (特殊) None 中性非极性基团，但归类于N，如某些胺的非极性部分 SN1a N S 1 Acceptor 弱中等极性氢键受体，如叔胺 (-N(CH3)-) TN4a N T 4 Acceptor 中等极性氢键受体，如醚氧 (-CH2†-O-CH2†-) Q0 Q R 0 (特殊) None 带形式电荷但无特定氢键能力的基团，或用于描述电荷离域的离子 SQd Q S (level varies) Donor 带电荷的氢键给体，如质子化的胺基 (-NH3+) TQa Q T (level varies) Acceptor 带电荷的氢键受体，如羧酸根 (-COO-) X3h X (R/S/T) 3 None h 含卤素化合物，如二氯乙烷中的 -CHCl-CH2Cl (X3h 代表一个氯原子和部分碳链) 注：上表仅为示例，具体的珠子类型选择和参数化应参考最新的 Martini 文档和相关研究。极性水平和氢键能力可能因具体的化学环境和参数化目标而有所调整。“N/A”表示不适用。 3. Martini 3 珠子参数化策略 (Martini 3 Bead Parametrization Strategy) Martini 3 珠子的参数化遵循系统性的方法，结合了“自上而下”（top-down，基于实验热力学数据）和“自下而上”（bottom-up，基于全原子模拟数据）的策略，旨在准确再现分子的物理化学性质 (2)。 3.1 非键相互作用 (Non-bonded Interactions) 非键相互作用的参数化主要目标是再现小分子在不同溶剂对之间的分配自由能 ($\Delta G_{\text{transfer}}$) (8)。常用的溶剂对包括正十六烷/水、正辛醇/水和氯仿/水等。通过拟合这些分配自由能，可以有效地校准溶质-溶剂以及溶剂-溶剂之间的交叉相互作用强度 (8)。第二个核心参数化目标是溶剂的互溶性数据，可以通过定性观察或计算混合过剩自由能来进行检验 (8)。互溶性数据同样能够反映不同种类分子间的交叉相互作用以及它们各自的自相互作用的相对强度。非键相互作用通常采用 Lennard-Jones (LJ) 势来描述： $V_{LJ}(r_{ij}) = 4 \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]$ 其中 rij 是珠子 i 和 j 之间的距离，$\sigma_{ij}$定义了珠子间的有效直径（相互作用为零的点），$\varepsilon_{ij}$定义了势阱深度，代表相互作用强度。对于带电荷的 Q 型和 D 型珠子，除了 LJ 相互作用外，还包含库仑相互作用： $V_C\left(r_{i j}\right)=\frac{1}{4 \pi \varepsilon_0 \varepsilon_r} \frac{q_i q_j}{r_{i j}}$ 其中$q_i$和$q_j$是珠子的电荷，$ε_0$是真空介电常数，$ε_r$是相对介电常数（在 Martini 中通常设为15，用于隐式地考虑水的屏蔽效应，但具体值可能因模拟体系而异）。 3.2 键合相互作用与几何中心映射 (Bonded Interactions and Center-of-Geometry (COG) Mapping) 与主要依赖实验数据的非键参数化不同，键合相互作用（键长、键角、二面角）的参数主要通过“自下而上”的方法获得，即参考全原子 (All-Atom, AA) 模拟得到的结构分布 (16)。 Martini 3 的一个重要改进是采用了基于“几何中心”(Center-of-Geometry, COG) 的映射规则来定义 CG 模型的键合参数，取代了 Martini 2 中常用的“质量中心”(Center of Mass, COM) 映射 (8)。COG 映射在计算分子片段中心时考虑了氢原子的位置，这使得 CG 模型能更好地保持其对应全原子参考结构的体积和形状 (8)。COM 映射有时会导致不满意的键长和过高的堆积密度，而 COG 映射则能产生更接近实际的分子性质（如溶剂可及表面积）和本体性质（如质量密度）(8)。对于接近全原子分辨率的映射（例如使用 T 型珠子），COG 映射尤为关键。例如，对于4对1映射的烷烃链，COM 和 COG 的差异不大；但对于2对1映射的苯环，两者差异显著 (8)。基于 COG 的键长可以直接从全原子模型中提取，这使得参数化过程更易于自动化。这些初始参数在需要更高精度时可以被进一步优化 (8)。常用的键合势函数包括：键长 (Bonds): 简谐势$V_b(l) = \frac{1}{2} K_b (l - l_0)^2$ 键角 (Angles): 简谐势$V_a(\theta) = \frac{1}{2} K_{\theta} (\theta - \theta_0)^2$ 二面角 (Dihedrals): 周期性势$V_d(\phi) = K_{\phi} [1 + \cos(n\phi - \phi_0)]$在某些情况下，特别是对于需要保持刚性平面结构或特定构象的分子，也会使用约束 (constraints) 或特殊势函数（如improper二面角）(8)。 4. Martini 3 映射方法学 (Martini 3 Mapping Methodology) 将全原子结构映射到粗粒化表示是构建 Martini 模型的首要步骤。Martini 3 提供了一套更一致的规则和指导原则，旨在优化 CG 模型的体积和形状表示。 4.1 基本原则 (Basic Principles) 进行原子到珠子的映射时，应遵循以下基本原则 (5)：原子分组: 通常将2-4个非氢重原子（及其相连的氢原子）映射为一个 CG 珠子。R、S、T 型珠子分别对应约4、3、2个重原子的映射。化学基团完整性: 尽量避免将特定的化学官能团（如酰胺基、羧基、完整的芳香环单元）分割到不同的珠子中，以保持其化学特性。对称性保留: 映射方案应尽可能尊重原始分子的对称性。体积与形状保持: CG 模型应能较好地再现全原子结构的体积和形状。COG 映射对此有重要贡献。珠子数量优化: 珠子的总数应被优化，目标是使每个珠子代表的重原子数与理想映射（如4:1, 3:1, 2:1）的最大偏差控制在每10个重原子中±1个非氢原子的范围内。原子共享: 在某些情况下，为了保持底层原子结构的对称性（例如在苯酚、四氢呋喃、甲苯等分子中），一个或多个原子可能被相邻的珠子共享（在一些文献的表格中用 † 标出）(8)。在从 COG 映射的全原子模拟中提取键合参数时，必须考虑到这种共享。 4.2 环状结构映射 (Mapping Ring Structures) 环状结构因其特殊的几何形状和电子特性，在 Martini 3 中有特定的映射策略，通常使用 S 型或 T 型珠子 (8)：芳香环 (Aromatic Rings): 芳香环（如苯环、萘环等）通常使用 T 型珠子进行描述，以更好地再现其平面性和堆叠行为。苯 (Benzene): 作为典型的芳香化合物，苯被模型化为三个 TC5 类型的 T 型珠子，每个珠子代表两个连续的碳原子及其相连的氢原子。TC5 是非取代芳香环中 -C=C- 基团的首选珠子类型。使用基于 COG 的键长（约 0.29 nm），可以很好地再现苯的液体密度。芳香环模型通常使用约束来连接珠子，因为其键长分布非常窄，需要非常刚性的势函数，这反过来又可能需要较短的模拟时间步长。对于更延展的刚性结构，可以考虑使用虚拟位点。脂肪环 (Aliphatic Rings): 脂肪环（如环己烷）通常使用 S 型珠子进行描述，以捕捉其相对于芳香环更大的体积感。环己烷 (Cyclohexane): 作为典型的脂肪环化合物，环己烷通常被描述为一个双 S 珠模型 (SC3-SC3)。SC3 珠子的选择基于分配数据。脂肪环模型通常使用键合相互作用（而非约束）连接，因为它们的键长分布相对较宽。使用约 0.378 nm 的键长，可以很好地再现环己烷的液体密度。 4.3 取代基与支链映射 (Mapping Substituents and Branched Chains) 对于带有取代基的环状结构或具有支链的分子，映射时需遵循以下两个主要原则 (8)：用最少数量的珠子映射所有非氢原子。尽可能保持分子的对称性、体积和形状，其中芳香环最好用 T 型珠子描述，脂肪环最好用 S 型珠子描述。例如：甲苯 (Toluene): 在苯环上增加一个甲基时，苯环原有的三个 T 型珠子中的一个会变成一个更大的 S 型珠子，以容纳额外的碳原子 (8)。乙苯 (Ethyl-benzene): 如果是乙基取代，则会为乙基额外增加一个 T 型珠子（代表乙基的两个碳原子），而苯环部分则可以保持其精确的三个 T 型珠子模型 (8)。支链烷烃/完全支化基团: 对于如新戊烷（包含5个非氢原子）这样的完全支化基团，通常会使用尺寸较小的珠子。例如，尽管有5个重原子，但由于中心碳原子被包埋，其对环境的暴露减少，因此可以使用一个 S 型珠子来代表整个新戊烷基团，或者根据具体情况进行更细致的划分 (8)。 4.4 高级模型设计策略 (Advanced Model Design Strategies) 对于具有多个芳香/脂肪环结构或复杂连接方式的小分子，Martini 3 提供了一些高级模型设计策略，常利用虚拟（相互作用）位点 (virtual sites) 来提高模型的数值稳定性和计算性能 (8)： “铰链”模型 (Hinge Model): 适用于刚性的稠合多环化合物，如萘 (Naphthalene)。一个简单的由5个 TC5 珠子通过约束连接的萘模型在凝聚相中可能导致数值不稳定。 “铰链”结构使用4个外部珠子，并将中心的一个或多个珠子描述为虚拟相互作用位点（其位置由构建粒子定义，受力会分配给构建粒子，质量均匀分配给构建粒子）。这种方法减少了约束数量，提高了数值稳定性和模拟速度。通常还会施加一个不当二面角来保持铰链模型的平面性。 “分而治之”模型 (Divide and Conquer): 适用于由刚性平面片段组成的任意长链，且需要控制片段间的相对二面角，这在小分子和共轭聚合物（如2,2’-联噻吩）中很常见。例如，两个噻吩环各用三个 T 型珠子描述。为了连接它们并控制二面角，可以在每个噻吩环的几何中心使用两个虚拟的非相互作用的哑位点 (dummy sites)，并通过简谐键连接这两个哑位点。然后可以在这些哑位点和每个噻吩环上的两个粒子（如硫原子）之间施加二面角势。 “分子转角”模型 (Molecular Turn): 用于处理通过 sp2 杂化碳连接的环系统，这种连接方式会产生一个“分子转角”（如某些具有特定扭转行为的分子）。为了保持扭转运动的正确旋转轴，需要特别注意。通常会使用虚拟哑位点：每个环体系的 COG 处各一个，第三个位于连接的 sp2 杂化碳上。通过在这些虚拟位点之间施加简谐键和角势，并辅以适当的二面角势和不当二面角势来维持正确的几何构型和旋转自由度。这些高级策略体现了 Martini 3 在处理复杂分子结构方面的灵活性和精确性。 5. 总结与拓展资源 (Conclusion and Further Resources) 5.1 总结 (Summary) Martini 3 力场通过引入新的珠子类型、更精细的尺寸划分 (R, S, T)、系统的命名规则（包含极性、氢键能力、电子特性等后缀）以及改进的参数化策略（特别是基于几何中心 COG 的映射），显著提升了粗粒化模拟的准确性和适用范围 (1)。其核心优势在于能够在保持较高计算效率的同时，捕捉到关键的化学物理特性，从而能够模拟更大尺度和更长时间尺度的分子过程。在对新分子（尤其是如本教程中所示的复杂聚合物）进行 Martini 3 映射时，关键步骤包括：仔细分析全原子化学结构，识别关键官能团。遵循原子分组（2-4个重原子/珠子）、化学基团完整性、对称性和体积/形状保持等基本映射原则。参考 Martini 3 珠子类型表和命名规则，为每个分子片段选择最合适的珠子类型和尺寸。对于缺乏直接预参数化珠子的特殊基团（如本例中的N-氧化物），需基于其化学物理特性（极性、氢键、分配行为等）类比选择最接近的现有珠子，或进行审慎的重新参数化。定义珠子间的键合连接。通过这些步骤，可以为目标分子构建合理的 Martini 3 粗粒化模型，为后续的分子动力学模拟打下坚实基础。 5.2 拓展资源 (Further Resources) 为了更深入地学习和应用 Martini 3 力场，以下资源非常宝贵： Martini 官方网站: http://cgmartini.nl (2)。这里可以找到最新的力场文件、教程、FAQ 以及已参数化的分子拓扑数据库（包括脂质、蛋白质、糖类、溶剂和小分子等 (8)）。主要出版物: Souza, P.C.T., Alessandri, R., Barnoud, J. et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. Nat Methods 18, 382–388 (2021). (1) (Martini 3 的奠基性论文)。 Alessandri, R., Souza, P.C.T., Thallmair, S. et al. A coarse-grained force field for small molecules: Martini 3. ChemRxiv (2021). (8) (针对小分子参数化的重要参考，包含大量构建模块信息)。模拟软件包: Martini 力场广泛应用于 GROMACS (4)。NAMD 等其他软件包也有相应的支持或转换工具 (26)。辅助工具: Martinize (或 Martinize2): 用于将全原子蛋白质（以及其他分子）结构转换为 Martini CG 模型的常用脚本 (5)。 Insane.py: 用于快速搭建复杂膜体系的脚本 (24)。 Polyply: 用于生成聚合物拓扑的工具 (29)。 MartiniGlass: 用于 VMD 中可视化 Martini 模型的 Python 包 (23)。力场参数下载: Martini 3 核心参数文件 (martini_v300.zip 或类似名称) 可从官方网站下载，其中包含了珠子定义 (.itp 文件)、相互作用矩阵以及多种已参数化分子的拓扑文件 (14)。小分子数据库通常托管在 GitHub 等代码仓库中，如 ricalessandri/Martini3-small-molecules (8)。 Marrink实验室的 GitHub 仓库 (marrink-lab/martini-forcefields) 也是获取最新参数和分子拓扑的重要来源 (29)。利用这些资源，研究者可以有效地将 Martini 3 应用于广泛的化学和生物物理问题研究中，探索复杂体系的结构、动态和热力学性质。

Molecular Dynamics · 2025-05-27

QM/MM Study of Enzymatic Reactions: Analysis of Chu et al. (2014) Paper

基本信息 Wen-Ting Chu, Qing-Chuan Zheng* and Hong-Xing Zhang 作者来自吉林大学理论化学研究所发表于 Phys.Chem.Chem.Phys., 2014, 16, 3946 DOI：https://doi.org/10.1039/C3CP53935K 论文摘要双磷酸甘油酸变位酶（bisphosphoglycerate mutase, BPGM）是一种多功能酶，其主要功能是合成血红蛋白的变构效应物——2,3-双磷酸甘油酸（2,3-BPG）。该酶亦可催化2,3-BPG水解生成3-磷酸甘油酸（3-PGA）。本研究通过量子力学/分子力学（QM/MM）方法，结合元动力学（metadynamics）和伞形采样（umbrella sampling）模拟，从理论角度揭示了人类双磷酸甘油酸变位酶（hBPGM）磷酸酶与合酶活性的反应机制。模拟结果不仅呈现了两类反应路径的自由能曲线，还阐明了活性位点中关键残基（如His11和Glu89）的作用。此外，反应能量势垒计算表明，hBPGM的合酶活性显著高于磷酸酶活性，且理论估算的势垒值与实验数据高度吻合。本研究为深入解析双磷酸甘油酸变位酶家族的催化机制提供了重要理论依据。关键词：双磷酸甘油酸变位酶；QM/MM模拟；自由能曲线；能量势垒；变构效应物 Introduction hBPGM是一种红细胞特异性多功能酶，具有合酶（EC 5.4.2.4）、变位酶（EC 5.4.2.1）和磷酸酶（EC 3.1.3.13）三种活性，其核心功能是催化1,3-双磷酸甘油酸（1,3-BPG）转化为2,3-双磷酸甘油酸（2,3-BPG）。作为血红蛋白的关键变构效应物，2,3-BPG通过稳定脱氧血红蛋白构象调控氧运输效率。尽管三种活性共享同一活性位点，实验表明合酶活性显著高于其他两种，而磷酸酶活性则负责水解2,3-BPG生成3-磷酸甘油酸（3-PGA）。Wang等人通过晶体结构研究（PDB: 2H4Z）揭示了活性位点残基His11与Glu89的催化作用，并提出磷酸酶反应遵循S2机制：His11作为亲核攻击位点夺取底物的磷酸基团，Glu89则通过质子转移稳定中间态。然而，hBPGM催化过程中原子尺度动态路径（如过渡态构型、自由能变化）仍缺乏理论解析。为此，本研究首次采用量子力学/分子力学（QM/MM）方法，结合元动力学（metadynamics）和伞形采样（umbrella sampling）模拟，系统性分析磷酸酶与合酶活性的反应路径与能量势垒，旨在从理论层面揭示hBPGM催化特异性的分子基础，为酶家族功能演化与药物设计提供新见解。 Fig. 1 The proposed mechanisms for the phosphatase and the synthase reactions. Methods 体系初始模型构建研究基于人源双磷酸甘油酸变位酶（hBPGM）与底物2,3-BPG的复合物晶体结构（PDB: 2H4Z，分辨率1.50 Å），选取单体链A（Ser2-Gln256）作为分子动力学（MD）模拟的初始结构。所有结晶水分子被保留，活性位点残基的质子化状态通过PROPKA在线工具（http://propka.ki.ku.dk/）确定：Glu89保持质子化，His11在δ位点单质子化，以匹配催化机制的需求。缺失的氢原子通过AMBER 12软件的LEaP模块添加，蛋白质参数采用ff99SB力场，底物2,3-BPG的参数由通用Amber力场（GAFF）生成。体系电荷通过添加钠离子中和，并置于TIP3P水分子填充的八面体周期箱中，确保蛋白质外层与水箱壁的最小距离为8.0 Å。分子动力学模拟 hBPGM/2,3-BPG复合物的经典MD模拟分为能量优化、平衡与生产三阶段：能量最小化：分两步进行，首先对水分子和离子进行2000步最速下降法+3000步共轭梯度法优化，随后对全体系重复相同流程以消除空间冲突。升温与平衡：在NVT系综下以1 K/ps速率升温至300 K，随后进行200 ps平衡模拟，期间对蛋白质Cα原子和配体原子施加弱限制（力常数0.5 kcal/mol/Å²）。 production模拟：在NPT系综下进行20 ns自由MD模拟，采用SHAKE算法约束氢键，粒子网格Ewald（PME）方法处理长程静电相互作用（截断值10 Å），时间步长2 fs。体系稳定性通过蛋白质骨架均方根偏差（RMSD≈1.3 Å）和配体构象（RMSD≈1.0 Å）验证，所有结构可视化由PyMOL完成。 QM/MM元动力学模拟基于平衡后的MD构象，采用AMBER软件结合PLUMED 1.3插件进行量子力学/分子力学（QM/MM）元动力学模拟。 QM区域包含底物2,3-BPG、His11和Glu89，MM区域为体系其余部分，QM/MM边界通过引入四个氢连接原子处理。每部分模拟运行1 ns，采用PM3半经验方法，高斯势宽度0.35 Å、权重0.1 kcal/mol，并设置±3.0 Å能量墙防止基团逃逸。自由能面（FES）通过累积的高斯势构建，过渡态（TS）通过能量最高点确定。在QM/MM元动力学模拟中，磷酸酶和合酶活性的反应路径通过原子间距离差作为集体变量（Collective Variables, CVs）进行描述，具体定义如下：磷酸酶活性第一步（磷酸基团转移至His11）反应坐标（ξ₁）：定义为底物磷酸基团的O3-P10键长与P10-His11的NE2原子键长之差，即： R(O3−P10)−R(P10−NE2) 物理意义：正值增大时，O3-P10键断裂（距离增大），P10-NE2键形成（距离缩短），反映磷酸基团从底物转移至His11的进程。第二步（Glu89质子转移）反应坐标（ξ₂）：定义为Glu89的OE2-HE2键长与HE2-O3（底物）键长之差，即： R(OE2−HE2)−R(HE2−O3) 物理意义：负值增大时，Glu89的HE2质子向底物O3转移，促进磷酸基团脱离（图1）。合酶活性第一步反应反应坐标（ξ₃）：定义为His11的P10-NE2键长与底物1,3-BPG的P10-O6键长之差，即： R(P10−NE2)−R(P10−O6) 物理意义：正值减小时，P10-O6键断裂（距离增大），P10-NE2键形成（距离缩短），反映磷酸基团从His11转移至底物的逆过程（与磷酸酶第一步相反）。后两步实际上就是磷酸酶催化的逆反应，不用再模拟一遍了。伞形采样验证为验证元动力学结果，对同一体系进行伞形采样分析。磷酸酶反应的两步及合酶反应的第一步被划分为多个窗口（步长0.1 Å，范围-3.0~3.0 Å），每个窗口进行50 ps采样（力常数200 kcal/mol/Å²）。初始构象从前一窗口末帧延续，采用PM3/ff99SB组合力场。数据通过加权直方分析法（WHAM）整合，去除谐波势影响后计算平均力势（PMF）。与元动力学相比，伞形采样在QM/MM边界处调整氢连接原子位置（Cα-Cβ键），以提高计算精度。 Results 普通MD模拟 hBPGM单体具有a/b折叠结构，包含两个域，六个β链和十个α螺旋。进行了20纳秒的MD模拟以获取该复合物的稳定构象，用于进一步机制研究。能量及稳定性评估总能量结果显示，在MD模拟后，复合物达到了平衡状态。蛋白质和配体相对于晶体结构的均方根偏差（RMSD）值表明，在整个MD运行过程中，蛋白质骨架RMSD稳定在约1.3 Å；而配体2,3-BPG在初始100皮秒后的RMSD保持在大约1.0 Å，没有发生构象变化。均方根波动（RMSF）分析显示蛋白质中有两个片段（Glu127到Gln151和Glu224到Gln251）存在较大的构象变化，但这些区域都是远离活性位点的柔性环区。氢键网络：2,3-BPG带五个负电荷并拥有十个氧原子作为氢键供体，与多个hBPGM残基形成了一系列氢键，包括Arg10、His11等。磷酸基团：2,3-BPG中的两个磷酸基团被不同的口袋包围，分别由特定的精氨酸和其他催化残基稳定，形成了反应中心，对于合成酶和磷酸酶活性至关重要。综上所述，通过MD模拟证明了hBPGM/2,3-BPG复合物已达到平衡，为后续的量子力学/分子力学（QM/MM）机制计算做好了准备。磷酸酶活性（Phosphatase Activity）的结果 hBPGM的磷酸酶活性催化2,3-双磷酸甘油酸（2,3-BPG）水解为3-磷酸甘油酸（3-PGA），其反应机制分为两步，通过量子力学/分子力学（QM/MM）结合元动力学（metadynamics）和伞形采样（umbrella sampling）方法进行模拟，具体结果如下： 1. 磷酸酶反应的两步机制与能量势垒第一步：磷酸基团转移（2,3-BPG → His11）反应坐标： ξ1监测O3-P10键断裂（距离从1.7 Å增至4.0 Å）和P10-NE2键形成（距离从4.0 Å缩短至1.8 Å）。能量势垒元动力学：25.75 kcal/mol（TS1b，对应ξ₁=0.81 Å）。伞形采样：21.61 kcal/mol（TS1a，ξ₁=-0.59 Å）。构象变化 His11的咪唑环旋转60°，形成共价键（图3B）。 Mulliken电荷显示O3电荷从-0.688（反应物R）变为-0.852（中间态I），NE2电荷从-0.178变为-0.102，表明电子重排（表1）。第二步：质子转移（Glu89 → O3）反应坐标： ξ2 监测Glu89的HE2质子转移至O3（OE2-HE2距离从1.2 Å增至2.3 Å，HE2-O3距离从3.0 Å缩短至1.3 Å）。能量势垒元动力学：5.21 kcal/mol（TS2，ξ₂=-0.1 Å）。伞形采样：6.32 kcal/mol（ξ₂=-0.18 Å）。 Glu89的作用 Glu89羧基旋转90°，与O3形成氢键，稳定中间态（图6）。 O3电荷从-0.930（中间态I）变为-0.434（产物P），OE2电荷从-0.353变为-0.701（表2）。 2. 方法比较能量势垒差异：伞形采样因更精细的窗口划分（步长0.1 Å）和氢连接原子优化（Cα-Cβ键），其势垒（21.61 kcal/mol）较元动力学（25.75 kcal/mol）更接近实验推算值（20.63 kcal/mol）。过渡态构象验证：两种方法的过渡态构型中，P10-NE2距离分别为2.4 Å（元动力学）和2.5 Å（伞形采样），高度一致（图5）。结论磷酸酶活性的限速步骤为第一步的高能量势垒（约20 kcal/mol），而Glu89的质子转移显著加速第二步反应。QM/MM模拟结果与Wang等人的实验数据（S2机制）一致，揭示了hBPGM催化中残基协同作用的分子基础，为靶向酶活性调控提供了理论依据。 His11：作为亲核攻击位点，直接参与磷酸基团转移。 Glu89：通过质子转移降低第二步势垒，促进磷酸基团脱离。 Arg10/Arg62：通过氢键稳定磷酸基团，降低反应能量需求（图2E）。合酶活性（Synthase Activity）的结果合酶活性催化1,3-双磷酸甘油酸（1,3-BPG）转化为2,3-双磷酸甘油酸（2,3-BPG），包含三个步骤，其中后两步为磷酸酶反应的逆过程。研究通过QM/MM元动力学和伞形采样模拟，揭示了以下关键结果： 1. 反应路径与能量势垒第一步（磷酸基团转移）：反应坐标定义为P10（磷酸基团磷原子）与His11的NE2原子距离差（ξ₃ = R(P10-NE2) – R(P10-O6)）。元动力学模拟显示能量势垒为12.98 kcal/mol（TS1），伞形采样结果为9.47 kcal/mol（图4B）。过渡态（TS1）对应ξ₃ ≈ -0.69 Å，此时P10-NE2距离从4.7 Å缩短至3.4 Å（元动力学）或2.5 Å（伞形采样），P10-O6距离从1.8 Å延长至4.5 Å（表3）。限速步骤（第三步：3-PGA → 2,3-BPG）：元动力学计算势垒为19.36 kcal/mol，伞形采样势垒为15.24 kcal/mol，与实验推算值16.49 kcal/mol（基于速率常数kcat = 13.63 s⁻¹）高度吻合（图7A）。产物态（2,3-BPG）自由能显著低于反应物态（-6.40 kcal/mol），表明反应热力学有利。 2. 原子相互作用与电荷变化 His11的动态作用 His11位于柔性loop区域，在第一步中向1,3-BPG移动并旋转约30°，捕获磷酸基团（图8）。 Mulliken电荷分析显示，NE2原子电荷从-0.253（反应物R）变为-0.126（中间态I1），O6原子电荷从-0.296变为-0.779，表明磷酸基团转移伴随电子重排（表3）。 Arg10与Arg62的稳定作用：这两个精氨酸通过氢键稳定磷酸基团，降低反应势垒。突变实验证实，Cys23和Ser24的突变（如C23T、S24G）显著降低合酶活性，因其破坏底物与蛋白质的氢键网络（图8）。 3. 构象变化与二面角调整底物构象重排：反应过程中，1,3-BPG的O5和O6原子向Cys23和Ser24旋转，形成新的氢键（图8）。二面角O5-C11-C1-O3从-27.93°（R态）变为82.07°（I1态），表明羟基（O3）向Glu89方向旋转，为后续质子转移做准备（表3）。后面就是O3被拔掉质子，夺回磷酸了 4. 方法比较元动力学 vs. 伞形采样：伞形采样因更精确的氢连接原子处理（Cα-Cβ键）和窗口划分（步长0.1 Å），其势垒值（9.47 kcal/mol）较元动力学（12.98 kcal/mol）更接近实验数据。合酶活性优势：合酶总势垒（15.24 kcal/mol）显著低于磷酸酶（21.61 kcal/mol），与实验测得的速率常数差异（合酶13.63 s⁻¹ vs. 磷酸酶0.0125 s⁻¹）一致，解释了hBPGM以合酶活性为主导的生理功能。评论：妙就妙在磷酸酶势垒最高的一步是N从O上抢走P，众所周知PO一家亲，而这正好为合酶提供了自由能的降低。人家合酶是拆掉磷酸-羧酸酐，自然势垒不那么高，还不用质子转移。结论合酶活性通过His11的定向移动、Arg10/Arg62的静电稳定及底物构象调整，高效催化磷酸基团转移。QM/MM模拟不仅验证了Wang等人提出的S2机制，还量化了残基协同作用对降低能量势垒的贡献，为设计调控2,3-BPG水平的药物提供了原子级理论依据。活性位点的其他残基图9展示了在磷酸酶反应的TS1b状态(A)和P状态(B)，以及合成酶反应的TS1(C)和R状态(D)下活性区域内的关键相互作用。与反应相关的残基被标记出来，它们之间的氢键以黑色虚线表示。这些信息强调了这些关键残基在催化过程中的重要作用。能量障碍：通常情况下，断裂一个O-P键需要大约80.06 kcal/mol的能量。然而，在hBPGM反应中，由于催化残基及其他活性位点残基的贡献，这一能量障碍显著降低。正电荷氨基酸的作用：底物磷酸基团周围存在多个正电荷的氨基酸（如精氨酸），它们与磷酸基团的负电氧原子有强烈的相互作用，有助于稳定过渡态。具体作用： Arg10 和 Arg62：两个精氨酸残基通过氢键与释放的磷酸基团相互作用，在磷酸酶和合成酶活性过程中帮助稳定过渡态。 His188：在反应开始前，His188通过氢键与His11相互作用，帮助其在hBPGM中达到正确位置并参与底物与蛋白质的结合。His188还形成氢键与磷酸基团及His11的NE2原子相连，减少了反应路径长度（从NE2到P10），使磷酸基团朝向His11的方向移动。在磷酸基团从底物转移到His11之后，Arg10、Arg62和His188通过氢键与磷酸基团相互作用，有助于保持中间体结构的稳定性。重要性：这些残基对于磷酸酶活性和合成酶活性反应至关重要，它们不仅降低了反应的能量障碍，而且通过特定的相互作用稳定了过渡态和中间体，从而促进了反应的进行。论文总结本文通过对人类bisphosphoglycerate mutase (hBPGM)的量子力学/分子力学 (QM/MM)模拟研究，成功地揭示了该酶在磷酸化和脱磷酸化反应中的催化机制，并提供了关于其动力学特性的定量估计。研究人员使用了经典分子动力学(MD)结合QM/MM和metadynamics以及umbrella sampling方法，这些方法为研究生物大分子的动力学行为提供了一种准确而有效的方法。研究结果表明，hBPGM的主要活性是合成酶，而不是磷酸酶或异构酶，这与之前的研究结果一致。通过本研究，可以更好地理解hBPGM在调节红细胞中2,3-BPG水平方面的作用，这对于深入研究hBPGM在疾病治疗方面的应用具有重要意义。未来展望未来可以通过进一步的研究，探索hBPGM在不同生理条件下的功能差异，例如氧气浓度、pH值等，以更好地了解其在调节2,3-BPG水平方面的作用。可以尝试将其他计算方法（如蒙特卡罗模拟）与QM/MM相结合，以更全面地研究生物大分子的动态性质。可以探索hBPGM与其他相关酶之间的相互作用，以更好地理解它们在代谢途径中的协同作用。个人Comments 比较经典的画反应过程的工作，也没有任何额外的东西，可以借鉴其流程。确定一个好的构象，开跑就完了表明His确实可以亲核进攻磷酸酯，能垒也确实挺高的，20多kcal/mol，看来我得be cautious了由于是用伞形采样画PMF，得到的是free energy surface；簇模型应该只能得到potential energy surface，看来还是MD好反正都是距离作为CV，伞形采样就行了，不用metaD了。metaD可能就是可以同时算多个距离，但伞形采样也可以设多个group吧。只设置距离作为CV，也能把二面角（单键转动）同时模拟出来部分结构图画得有点丑。。。不能把蛋白残基和底物区分一下颜色嘛看来画FES就应该是直接用WHAM算出来的结果作图，没想象中这么麻烦文章内容总结主要由AI完成，如有错误恳请指出！

Molecular Dynamics · 2025-03-13

File Conversion Among MD Simulation Engines Using ParmEd

File Conversion Among MD Simulation Engines Using ParmEd ParmEd is a versatile Python library that facilitates the interconversion of files between popular molecular dynamics (MD) simulation engines like Gromacs, Amber, and NAMD (CHARMM). This tool is especially useful for researchers and students working in molecular dynamics who need to switch between simulation packages without hassle. For example, you want to avoid setting up a protein-ligand complex in Gromacs (adding ligands to gmx force field files can be troublesome!) but do want to run MD simulations in Gromacs for its speed. You will need to use ParmEd to convert the Amber files to Gromacs format. Note that the MD engine uses different algorithms and settings. You cannot either adopt special settings in another MD engine (e.g. restraints, you should set it up again). You should not even wish to fully replicate a Gromacs simulation in Amber. But for most biological systems (e.g. the solvent is not that important), MD engine usually affects your simulation much less than other options, like the choice of force field. So feel free to switch between MD engines! Jump to the code section if you want a solution only. Installing ParmEd Here’s how you can install ParmEd using Anaconda: conda install -c conda-forge parmed If you have compiled Amber on your system, you might already have ParmEd installed as part of the AmberTools suite. To ensure it is properly integrated, refer to the comprehensive guide on compiling Amber, which is particularly useful if you are setting up everything from scratch. Introduction Knowing the file formats These file formats are what we need in MD simulations: Engine Construction Tool Topology file Coordinate file Parameter file Gromacs pdb2gmx .top/.itp .gro – Amber tleap .prmtop .inpcrd – NAMD VMD psfgen .psf .pdb .prm ParmEd logics ParmEd works simply: read in the topology and coordinate files, and write out two files in the desired format. ParmEd writes the parameters into .inpcrd (as it is) and .top files. Always find .prm files when converting both from and to NAMD. Other You can edit the system in ParmEd, which is out of the scope of this post. The file parsing is very detailed so you can manipulate the system as you like. Consult the ParmEd documentation for more details. Conversion Code The following code shows a framework of file conversion. It implements the basic residue renumbering function: you can set the starting residue number. The command is python xxx.py <system_name> <starting_residue_number> Your topolgy and coordinate files should be named <system_name>.xxx both. Note that we use offset-1 in the code since by default ParmEd residue numbers start from 1. [!WARNING] Always double check after the conversion! [!WARNING] For a very large system (hundreds of thousands of atoms), this process could take some time. From Amber to Gromacs # python amber2gmx_via_parmed.py pro 689 import parmed as pmd import sys prefix = sys.argv[1] offset = int(sys.argv[2]) amber = pmd.load_file(prefix+'.prmtop', prefix+'.inpcrd') # renumbering for residue in amber.residues: _ = residue.idx # Get the original index residue._idx += offset-1 residue.number += offset-1 # Save the modified files in Gromacs format amber.save(prefix+'.top', overwrite=True, combine='all') amber.save(prefix+'.gro', overwrite=True, combine='all') Gromacs sub-topology .itp files can be read, but cannot be written, i.e. ParmEd writes huge topology/coordinate files without subfiles as in Amber/NAMD. From CHARMM to Gromacs # python charmm2gmx_via_parmed.py pro 689 import parmed as pmd from parmed.charmm import CharmmParameterSet import sys prefix = sys.argv[1] offset = int(sys.argv[2]) structure = pmd.load_file(prefix+'.psf') # renumbering for residue in structure.residues: _ = residue.idx residue._idx += offset-1 residue.number += offset-1 parameter = CharmmParameterSet('par_all36m_prot.prm', 'toppar_water_ions_namd.str') # add more if necessary # edit the sign of epsilon for atomname, atomtype in parameter.atom_types.items(): atomtype.epsilon *= -1 atomtype.epsilon_14 *= -1 structure.load_parameters(parameter) # Save the modified files in Gromacs format structure.save(prefix+'.top', overwrite=True, combine='all') structure = pmd.load_file(prefix+'.pdb') structure.save(prefix+'.gro', overwrite=True, combine='all') [!TIP] ParmEd does not realize that for epsilon gmx adopts the absolute value while charmm files store the real value (negative!) [!NOTE] In parameter files like par_all36m_prot.prm downloaded from CHARMM website, officially all atom type definitions are commented, but we should uncomment them for parmed, or it cannot find atomtypes. Or read .rtf files too. Double check your files! From Gromacs to Amber # python gmx2amber.py system import parmed as pmd import sys prefix = sys.argv[1] parm = pmd.load_file(prefix+'.top', prefix+'.gro') # Save the modified files parm.write(prefix+'.prmtop') parm.write(prefix+'.inpcrd') I actually have not tried this (see problems). You may need to add residue renumbering mechanisms. Practice yourself! And I guess from CHARMM to Gromacs works similarly. Renumber gmx files This adopts the similar process. The original files are overwritten. # python gmx_renumber_via_parmed.py pro 689 import parmed as pmd import sys prefix = sys.argv[1] offset = int(sys.argv[2]) gmx = pmd.load_file(prefix+'.top', prefix+'.gro') # renumbering for residue in gmx.residues: _ = residue.idx residue._idx += offset-1 residue.number += offset-1 # regenerate and revalidate the internal parameters, usually do this after modifying the structure gmx.remake_parm() # Save the modified files gmx.save(prefix+'.top', overwrite=True) gmx.save(prefix+'.gro', overwrite=True) From CHARMM to Amber To convert CHARMM files to Amber format, use chamber: chamber -top topol.rtf -param params.par -str stream.str -psf structure.psf -crd structure.crd -outparm amber.prmtop -outcrd amber.inpcrd Topology files (-top, -str) are only necessary if the parameter files do not define the atom types Parameters (-str, -param) are applied to your structure -crd option accepts file formats like PDB, CHARMM CRD, Amber restart, etc. Issues Residue renumbering Problem None of these file formats are perfect. Gromacs files do not have chain identifiers. By default chains are separated into a few .itp files, so it’s hard to locate an atom in a specific chain in a .gro file. Amber files always start with residue numbers 1, which causes trouble when aligning with the “biological” residue nubmers. VMD files have full identifiers. However, we have to manually separate the chains when modeling. You cannot change the file formats unless your write your own MD engine. So just put up with it… With ParmEd, you can try to edit the residue numbers to match the “biological” residue numbers. Sadly, if you have multiple chains and they are overlapping, you still have to use that sequential residue numbers. But if you have only one chain, this won’t bother you. Edit in VMD During visualization in VMD, you can edit the residue numbers like this: mol new system.prmtop type parm7 first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all mol addfile md.nc type netcdf first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all # select whatever you are interested, but too many water many slow down the process set all [atomselect top "protein or resname LIG or resid 1 to 1500"] foreach idx [$all get index] { set atom [atomselect top "index $idx"] $atom set resid [expr [$atom get resid] + 688] } Edit in ParmEd In ParmEd, every Residue object in a Structure has an idx attribute. This attribute indicates the residue’s index within the structure, and it is managed internally by ParmEd. It is crucial not to modify this attribute directly, as it could lead to inconsistent state within the structure. Some other attributes are also private and cannot be modified. Anyway, I’ve figured out the code to edit residue numbers. I don’t really know why I have to manipulate _idx, but it works. Feel free to inspect the attributes when debugging in your IDE, and create your own workflow! Parameters and atomtypes GROMACS: Independent Parameter Specification In GROMACS, topology files (typically .top) allow for each bond term to be specified independently. This means that different bond parameters can be assigned to the same pair of atom types, provided they occur in different contexts within the molecule. Example of a GROMACS bond specification: ; Bond parameters ; i j func length force_const 1 2 1 0.123 456.7 ; Asymmetric bond A 2 3 1 0.123 456.7 ; Asymmetric bond B CHARMM: Type-Based Parameter Definition Conversely, CHARMM typically defines parameters between different atom types based on a consistent set of parameters across all bonds involving those atom types. This approach assumes that identical pairs of atom types will always exhibit the same bonding characteristics, regardless of their molecular environment. BONDS CA CB 340.0 1.529 ; Standard peptide bond CA CG 317.0 1.510 ; Standard alkane bond Resolving Parameter Inconsistencies When converting from GROMACS to CHARMM formats using tools like ParmEd, discrepancies in how bond parameters are specified can lead to errors. For instance, ParmEd might encounter a ParameterError if it detects different bond parameters for the same atom types, which is permissible in GROMACS but not in CHARMM. This issue is particularly evident with complex ions or molecules optimized asymmetrically through QM methods, such as $\ce{Al(OH)(H2O)5^2+}$. To address these conversion challenges, users have two main options: Assign Different Atom Types: Modify the topology to assign unique atom types for bonds that require different parameters. Uniform Bond Parameters: Standardize bond parameters for each pair of atom types, ensuring consistency across the entire molecule. For more details on handling these conversions and the underlying code structure of ParmEd, consider exploring the following resources: ParmEd GitHub repository Issue related to parameter mismatches Discussion on handling different parameters End We welcome your feedback and contributions! If you have developed new workflows or if you encounter any issues, please don’t hesitate to reach out. For reporting problems, consider opening an issue on the ParmEd GitHub repository. Your insights and experiences are invaluable in enhancing the tools and community resources.

Molecular Dynamics

Mendelevium

Contact

Molecular Dynamics